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DCTDOS : NEUTRON  AND  GAMMA  PENETRATION  IN  COMPOSITE  DUCT  SYSTEMS 


L.  V.  Spencer 

National  Bureau  of  Standards 
Washington,  D.C.  20234 


I . Introducti on 

This  paper  describes  computer  methods  for  estimating  neutron  and 
gamma  ray  fluence  rate,  dose,  and  even  spectral  features  due  to 
penetration  through  a series  of  duct  segments.  The  procedure  links 
together  data  for  individual  segments  --  straight  sections  and  bends  -- 
in  arbitrary  combinations;  and  the  resulting  composite  can  include 
computations  for  a room  at  the  end,  if  there  is  one.  This  particular 
method  was  developed  for  rapid  estimates  in  the  protection  problems 
against  nuclear  weapons;  but  the  concepts  which  are  employed  are  more 
broadly  appl  icabl  e. 

The  study  of  radiation  penetration  in  ducts  and  other  cavities  has 
developed  a large  literature,  which  can  be  subdivided  into  a)  studies 
of  neutron  and  gamma  ray  albedos,5’13’18-24’29’38’43-45’58’61’78,  b) 
experiments, l»3-4»6»8-9»12»15-16»27-28»30»33-34»37»39-40»48»54’59» 

70-71,  c)  Monte  Carlo  duct  analysis,17’25’33-37’43’46"50’52"59’75,  and 
d)  some  form  of  analytical  analysis.7’12"13’26’31-32’41’53’60’62-70 
Summary  discussions  of  the  various  problems,  data,  and  methods  are  to 
be  found  in  many  places,  among  which  are  Refs.  15,41-42,55-58,  and  74. 

Many  reactor  shielding  studies  are  not  referenced  above. 


The  program  described  here  concludes  a series  of  studies  which 
started  in  the  middle  and  latter  1960's.  In  these  studies  various 
aspects  of  radiation  penetration  through  ducts  were  examined  in  regards 
to  the  possibility  of  viewing  complex  duct  configurations  as  composites 
of  elementary  segments.67"69  Many  ideas  appearing  in  this  paper  were 
described  previously  in  Ref.  69.  While  this  approach  fits  mainly  with 
category  d)  above,  it  is  neither  a quadrature  nor  a curve-fitting 
method,  nor  does  it  emphasize  center-line  dose. 

This  report  presents  a usable  computer  program  (DCTDOS)  based  on 
concepts  described  in  that  earlier  study  and  gives  applications  to  some 
well  known  experimental  duct  studies.  Several  of  our  basic  assumptions 
are  as  follows: 

■k 

1)  That  the  relevant  flux,  which  depends  on  both  space  and 
direction  variables,  can  be  usefully  described  by  introducing  the 
order-of-scatteri ng  index  n (or  N=n+1).  In  DCTDOS  we  include  up  to 
n = 20,  which  is  barely  large  enough  for  complex  systems  with  large 
reduction,  and  may  eventually  be  enlarged. 

2)  The  albedo  angular  distribution  for  neutrons  is  well 
approximated  by  the  elementary  distributions  of  Chandrasekhar 1 1 for 
isotropic  scatterers. 69  Since  the  sensitivity  of  most  problems  to 
details  of  this  angular  distribution  is  rather  low,  we  use  it  to  obtain 
solutions  to  elementary  spatial  configurations  in  terms  of  the  variable 
ru  That  is,  we  use  one-velocity,  i sotropic-scatterer  albedos  to 
obtain  data  for  simple  segments  as  functions  of  ru 

3)  The  neutron  flux  component  which  has  been  backscattered  from 
concrete  walls  (say)  a total  of  n times  has  a spectral  distribution 
(for  a given  radiation  source)  which  is  assumed  to  depend  on  energy 
only  through  _n.  This  takes  advantage  of  a near-factori ng  of  the  single- 
reflection albedo  into  a directional  part  and  an  energy  part.  67-69  But 

Here  and  throughout  this  manuscript  the  term  flux  stands  for  fluence 
rate.  (See  ASTM  Standard  E170,  1986  Annual  Book  of  ASTM  Standards, 

Section  12.)  Flux  can  also  designate  a time  integral,  i.e.,  a fluence 
if  the  source  is  likewise  integrated  over  time. 
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in  turn  this  leads  us  to  assume  a factoring  of  the  whole  problem  into 
space-and-di rection  studies,  and  mul ti -refl  ection  energy  studies. 

Tables  of  mul ti -refl  ection  energy  albedos  as  functions  of  both  initial 
and  eventual  spectral  energy  are  required.69 

4)  In  using  the  output  of  one  segment  as  input  into  another,  some 
information  about  the  direction  of  particle  currents  is  a requirement. 
We  have  reduced  this  requirement  to  its  least  restrictive  form  by  using 
only  a single  parameter:  For  this  we  use  the  mean  reciprocal  cosine  of 

the  obliquity,  and  which  we  refer  to  in  the  convenient  form 
(-1  + 2/<sece>)  = s.  Note  that  <sece>  is  also  the  ratio 
fl ux/current . We  record  and  use  s as  a matching  parameter, 
distinguishing  between  unreflected  and  reflected  components  in  the 
preceding  of  any  two  segments  being  matched,  since  these  can  be  very 
di fferent . 

(5)  Radiation  reflected  in  cylinders  of  circular  cross  section, 
spheres,  and  between  infinite  plane  surfaces  has  the  property  of 
reentrance  at  the  emergence  obliquity  when  travel  in  the  wall  can  be 
ignored.  Hence  mul ti -refl ecti on  albedos  in  these  configurations  are 
identical  if  the  initial  incidence  distribution  and  spectrum  are 
identical.  We  take  advantage  of  the  fact  that  cylinders  with  circular 
and  square  cross  sections  have  a rough  equivalence,  as  do  cubes  and 
spheres.  Thus  we  believe  that  in  a fairly  usable  first  approximation 
mul ti -refl ection  distributions  between  infinite  plane  surfaces  has  a 
general  applicability  to  far  more  complex  configurations. 

Combination  of  segments  proceeds  by  simple  convolutions  of 
successive  functions  of  n.  Spectra  and  different  types  of  dose  are 
extracted  by  additional  convolutions  with  mul ti -refl  ecti on  energy  data. 
Hence  DCTDOS  marches  to  its  solutions  by  a series  of  convolutions  over 
data  interpolated  for  elementary  segments.  As  may  be  readily  guessed, 
such  a process  is  fast  and  can  be  flexible. 

Most  of  the  data  on  duct  sections  which  DCTDOS  utilizes  were 
developed  in  a series  of  unpublished  studies  by  S.  Woolf.76”77  The 
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mill  ti -refl  ection  albedo  energy  data  for  neutrons  is  due  to  Hagan  and 
Simmons,38  and  correspondi ng  data  for  gamma  rays  has  been  computed  by 
Spencer  and  Woolf.78  Since  all  data  are  rather  simple  in  concept  and 
readily  obtained  on  file,  we  anticipate  supplementation  and  replacement 
if  these  methods  are  further  developed. 

Gamma  ray  albedos  are  more  complex  than  neutron  albedos,  with  strong 
dependence  on  energy  above  0.5  MeV  [and  for  single  scattering].  If  the 
single-scattering  is  given  special  attention,  there  is  merit  in  treating 
gamma  rays  as  if  they  were  neutrons,  but  using  mul  ti -refl  ecti on  albedo 
energy  data  for  gamma  rays  specifically.  This  simple  approach  turns  out 
to  work  rather  well.  This  confirms  that  the  n = 1 reflection  component 
for  gamma  rays  dominates  the  difference  from  the  correspondi ng  neutron 
angular  distribution,  and  also  that  the  grazing  incidence  case, 
sometimes  termed  "skipping,"  is  comparatively  unimportant.  Differences 
of  albedo  magnitude  in  the  basic  segment  data  can  be  roughly  corrected 
by  use  of  simple  ratios  raised  to  the  n'th  power.  Recomputation  of  the 
space-angular  data  for  gamma  rays  specifically  has  not  yet  been  done. 

In  Section  V we  present  results  for  both  neutrons  and  gamma  rays  in 
comparison  with  experimental  data  and  results  from  well  known  and 
detailed  computations:  for  neutron  data  in  the  Maerker  and 

Muckenthaler  duct46-49  and  for  gamma  ray  data  in  one  of  the  AMF 
ducts.73  We  discuss  spectral  and  order-of- refl ecti on  comparisons  as 
well.  A high  degree  of  precision  is  neither  attained  nor  expected;  but 
we  feel  that  absolute  values  and  trends  are  sati sfactori ly  reproduced 
for  both  types  of  radiation. 

We  have  organized  this  report  to  put  practical  information  necessary 
to  apply  DCTDOS  towards  the  front,  with  explanatory  and 
supplementary  material  located  towards  the  rear.  The  comparison  of 
results  for  classic  experiments  constitutes  a transition  section 
between  these  two  types  of  information. 
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II.  Structure  of  DCTDOS 


In  the  diagram  of  Fig.  1,  the  capitalized  words  represent 
subroutines,  the  function  of  each  being  suggested  by  choice  of  name. 

Data  input  are  of  three  different  types.  (INPUT)  represents  problem 
parameter  data.  The  basic  data  on  which  the  program  depends  are 
introduced  by  calling  ENTDAT  to  set  up  elementary  segment  data  and 
SRSDET  for  spectra  and  mul ti -refl ection  response  weights,  both  making 
use  of  types  of  mul ti -refl ection  albedo  energy  data. 

Parameter  input  begins  with  new  Run  data,  and  later  is  indicated  by 
four  parallel  lines  designated  Run,  Problem,  Length,  and  Weights.  A 
new  Length  changes  only  a single  (configuration)  parameter,  while  a new 
Run  brings  in  all  the  data  necessary  for  a new  duct  computation.  A new 
Problem  uses  previous  source  and  detector  information,  while  "New 
Weights"  changes  source  and/or  detector  while  retaining  the  previous 
confi guration. 

Branch  points  have  been  identified  by  lower-case  letters  a_  through  i_, 
their  significance  being  as  follows: 

a)  Choice  between  DCTISO  (isotropic  source  flux)  and  DCTSLT 
(mono-obliquity  source  angular  distribution). 

b,c,d)  The  option  here  has  to  do  with  completion  vs.  an  additional 
segment  of  duct,  according  to  the  complexity  of  the 
confi guration. 

e)  The  configuration  here  may  call  for  room  computations  or  for 
dose  computations  in  the  duct. 

f)  Choice  of  mul  ti -refl  ecti  on  albedo  energy  data  and  application 

of  these  data  to  the  dose  weights:  neutrons  (SND)  or  gamma 

rays  (SGD). 

g)  This  option  permits  recomputation  with  different  dose  weights. 
The  rerun  through  SRSDET  gives  wide  latitude  to  this  option. 

h)  This  option  permits  modification  of  one  of  the  straight  segment 
lengths,  followed  by  recalculation  with  all  other  parameters 
fixed . 

i)  The  choices  here  are  for  a new  configuration  with  unchanged 
source  and  detector  data  (16)  or  a fully  new  situation  (New 
Run,  13).  Termination  is  also  an  option  at  this  point. 
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III.  Input  Data  and  Options 


We  use  the  term  "card"  here  with  reference  to  a line  of  input  data, 
which  formerly  was  punched  on  an  IBM  card. 

In  preparing  for  a DCTDOS  run,  several  different  types  of  data  card 
are  required.  The  first  specifies  mainly  source  and  detector  data  to 
be  used,  and  we  refer  to  this  as  a Run-Description  card.  Following 
this  is  a card  with  index  numbers  mainly  for  specifying  outputs,  new  or 
old  dose  weights,  and  a series  of  segment  lengths  or  detector  positions, 
together  with  a title  card.  We  refer  to  these  as  Probl em-Descri pti on 
cards.  Finally,  there  is  a set  of  cards  giving  dimensions  for  a 
sequence  of  duct  segments,  which  we  refer  to  here  as  Conf i gurati on 
cards.  Sample  run  decks  are  indicated  in  Figs.  2 and  3;  here  we  mainly 
give  descriptions,  together  with  the  options  currently  available. 

A.  Run-Description  Card.  The  first  run  card  specifies  five 
integers,  followed  by  a decimal  number  in  exponential  form,  using  the 
Fortran  format  ( 515, 1 PE10 .3) . The  numbers,  in  their  order,  have  the 
following  interpretations: 

ISORS  specifies  the  type  of  radiation  spectrum  selected, 

IDOSE  specifies  the  type  of  detector  response  selected, 

ING  specifies  gamma  rays  or  neutrons,  with  the  inclusion  or 
dismissal  of  wall -capture  gamma  rays, 

LPIPE  selects  from  several  types  of  (straight)  duct  data  available 
for  neutrons  (only  a single  type  of  gamma  ray  data  is  used), 
IEL  first  specifies  the  type  of  output  by  quick  identification 
with  IPRINT.  Later  this  parameter  has  a variety  of  other 
uses  through  its  COMMON  storage  in  most  subprograms. 

DIMEN  is  a numerical  factor  which  can  change  or  rescale  output 
units. 
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Table  1 lists  options  currently  available  according  to  values  of  the 
five  index  numbers  referred  to  above.  But  special  spectra  or  response 
functions  can  be  readily  inserted,  either  by  adding  to  the  lists  above, 
or  by  insertion  as  a replacement  of  one  of  the  cases  listed  in  Table  1. 
We  have  found  this  useful  in  substituting  137Cs  source  data,  say,  for  a 
60Co  gamma  spectrum. 


Table  1.  Options  currently  available,  as  specified  by  the  different 
input  integers. 

ISORS 

1:  Neutron  INR  spectrum  (neutrons/group/incident  dose),  see  Ref.  2, 
Table  10) 

2:  Thermal  neutron  spectrum  (normalized  to  unity) 

3:  Epithermal  neutron  spectrum  (neutrons/sec/watt/cm2)  incident.  See 
Ref.  48,  Fig.  2,  based  on  34,200  n/sec/watt  into  (91.44  cm)2. 

4.  Gamma  ray  INR  spectrum  (gammas/group/incident  dose).  See  Ref.  2, 
Table  10). 

5.  60Cobalt  gamma  rays  (normalized  to  unity).  50  Curies,  2 photons 
per  interaction;  but  only  half  moves  forward  into  a duct  of 
6'x6'cross  sectional  area.  This  gives  3.319E9  gammas/cm2/mi n . 

IDOSE 

1:  Snyder-Neufel d dose  ( rads/( n/cm2) ) 

2:  Unit  weights  (for  flux  computations,  e.g.  neutrons/ ( n/cm2) ) 

3:  Henderson  tissue  dose  data  ( rads/ ( n/cm2) ) 

4:  Thermal  neutron  response  only  (neutrcns/( n/cm2) ) 

5:  Gamma-ray  absorbed  dose  data  (rads/(gamma/cm2) ) 

ING  LPIPE 

-1:  Wall -capture  gammas  only  1:  Finite-duct  (exit)  neutron  data 

0:  Neutrons  only,  or  gammas  only  2:  Semi -i nfi ni te-duct  data,  forward 
1:  Neutrons  + wall -capture  gammas  3:  Ditto,  both  forward  and  return 

1-3:  Any  value  for  gamma  rays 
IEL 

0:  ( I PR  I NT ) - - No  print-out,  first  case; 

1-4:  ( I PR I NT ) — Print-out  (Table  2,  also  see  discussion  in  Section  IV) 
DIMEN 

This  constant  factor  can  change  dimensions  for  the  weights  and/or 
renormalize  the  input  current. 


B.  Problem  Description  Cards.  Let  us  first  dismiss  the  title  card, 
which  is  the  last  of  these  cards.  Its  format  is  presently  (55H  ...  ). 
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The  first  input  card  has  again  the  format  ( 515 , 1PE10 .3) . The 
integers  specified  have  the  options  and  i nterpretations  listed  in 
Table  2 below,  with  entries  in  the  same  order  as  on  the  card. 

Following  the  main  input  card  is  a card  with  a single  integer 
specifying  the  number  (JDOS)  of  alternative  sets  of  weights.  And  this 
is  followed  by  JDOS  alternative-weights  cards,  each  again  using  the 
format  (515, 1 PE10 . 3 ) . These  cards  list  values  for  IPRBLM,  ISORS, 

IDOSE,  ING,  LP,  and  DIMEN  in  that  order.  It  has  turned  out  that  of 
these,  IPRBLM  is  a dummy  without  application,  and  can  be  equated  to  its 


value  on 

the  first  input  card. 

Table  2. 

Options  for  the  "Problem"  Parameters 

IPRBLM: 

-1,-2,- 

1:  The  after  this  after  this  one  is  to  use  the  same  run- 
weights  card. 

2:  The  computation  after  this  one  will  have  a new  run-weights 
card. 

3:  This  is  the  last  computation. 

-3:  The  present  computation  will  use  weights  already  computed 
and  stored  from  the  run-weights  card  which  precedes. 

NEL: 

This  gives  the  number  of  duct  elements  (plus  one,  for  the 
source).  This  is  also  the  number  of  case  cards  which  follow. 

IPRINT: 

Acceptable  values  are  0,  1,  2,  3,  4.  Zero  indicates 
no  printout;  1 means  printout  of  input  information 
together  with  computed  dose  and  dose  ratio;  2 also  gives 
weights  plus  dose  information  at  earlier  junctions;  3 adds 
group  spectra  and  mean  secant  directional  information;  and  4 
adds  intermediate  printout  from  the  subroutines  for  use  in 
checking  parameter  sizes  and  other  inquiry. 

JVARY: 

This  gives  the  index  ( IC)  for  the  straight  section  whose 
length  variation  may  be  desired.  Thus,  JVARY  = 2 would  refer 
to  the  first  straight  section  of  a duct. 

JDEL: 

This  index  gives  the  number  of  values  which  the  straight 
section  specified  by  JVARY  is  to  have.  JDEL  = 1 means  that 
only  one  value  is  to  be  computed,  the  most  common  case:  The 

values  of  JVARY  and  DEL  are  then  not  relevant  and  DEL  may  be 
omitted  or  set  to  zero. 

DEL: 

This  fixed  point  number  gives  the  increment  by  which  the 
length  of  the  JVARY  segment  is  to  be  increased  for  each  of 
JDEL  cases.  It  is  important  that  DEL  is  length, 
i.e.,  length  divided  by  the  square  root  of  the  cross  sectional 
area. 
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C.  Configuration  Cards.  These  cards  in  sequence  describe  the 
configuration.  We  discuss  them  in  terms  of  the  names  of  the  lists  into 
which  they  are  put  (in  DATDCT).  These  lists  are  NEXT,  WDTH,  HGHT, 

LNTH,  OTHR,  SELECT:  The  Format  used  is  (15,  5F10.3)). 

The  first  entry,  NEXT,  gives  IEL  values  for  a series  of  subprograms 
to  be  called  in  the  sequence  given  (DCTDOS,  order  35).  These  are 
listed  in  Table  3 below  for  easy  reference: 


Table  3.  NEXT  indices  used  to  call  the  main  subprograms. 
IEL  SUBROUTINE  Type  of  Section 


0 

Descriptive  parameters,  only  used  with  IC  = 1 

1 

IXTISO 

Straight  segment,  isotropic  source 

2 

DCTSLT 

Straight  segment,  slant  source 

3 

TURN 

Bend  segment 

4 

BTODR 

Straight  segment  occurring  after  a bend 

5 

BTODR 

Straight  segment  after  a bend 

6 

R00MX 

Room  segment 

7 

DOSE 

Applies  weights  and  computes  total,  detector  in  or  at 
the  end  of,  a bend  or  straight  section 

8 

OOSX 

Applies  weights  and  computes  total,  detector  pattern 
in  a room 

The  remaining  entries  on  the  case  description  cards  include  widths, 
heights,  and  lengths.  We  summarize  these  in  Table  4 using  the 
abbreviations,  W,  H,  and  L,  respecti vely . We  use  dimensioned  numbers, 
usually  in  feet,  in  practice;  the  program  does  its  own  scaling.  Those 
parameters  which  are  not  W,  H,  or  L require  special  explanations,  which 
fol 1 ow  Tabl e 4. 

Table  4.  Input  data  assignments  for  the  case  description  cards.  (Note 
Appendix  B.)  NEXT  corresponds  to  IEL  in  Table  5. 


NEXT 

WDTH 

HGHT 

XLTH 

OTHR 

ELECT 

I ZD 

SO 

PEN 

ZD00R 

RMCASE 

1 

W 

H 

L 

ARAT 

S = 0 

2 

w 

H 

L 

ARAT 

S = -1,0,1 

3 

w 

H 

RCSS 

F 

S = -1,0,1 

4,5 

w 

H 

L 

ARAT 

S = -1,0,1 

6 

W,WTH 

H 

L 

WD 

S = -1,0,1 

7 

RDLNTH 

S = -1,0,1 

8 

WDP 

NW-1 

NL-1 

AEX 

S = -1,0,1 
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The  following  interpretations  apply  to  the  various  quantities  in 
Tables  3 and  4. 

SELECT(IC)  is  normally  zero  or  omitted.  If  it  equals  1.,  the  reflected 
component  is  equated  to  zero  for  computations  by  the 


subroutine  selected  by  NEXT ( IC ) . If  it  equals  -1.,  the 
unreflected  component  is  equated  to  zero.  If  0.,  the  normal 
case,  both  contribute.  Note  that  "unreflected"  and 
"reflected"  here  refer  only  to  the  preceding  section  and  not 
to  N = 1 and  N > 1 in  general. 

SO 

This  index  is  zero  for  the  cosine  (current)  source.  In 
general  this  quantity  is  = -1  + 2./<sece0>  for  slant 
obliquity  90. 

PEN 

is  an  assumed  effective  scattering  depth  in  side-walls  prior 
to  reemergence.  It  has  normally  been  set  equal  to  zero;  its 
dimensions  are  those  of  W,  H,  and  L. 

ZDOOR 

is  a door  or  wall  thickness  assumed  (for  IZD>0)  at  the 
entrance  to  the  shelter.  Note:  It  is  limited  to  dimensions 

of  feet  of  concrete  presently,  and  assumes  no  modification 
of  source  spectrum  or  directional  distribution. 

RMCASE 

has  values  +1.,  0.,  -1.:  The  value  +1.  is  used  if  the  flux 

entering  the  room  (i.e.,  R00MX)  and  unreflected  by  the 
preceding  segment,  is  directed  towards  a nearby  sidewall. 

The  value  0.  is  assigned  if  this  flux  component  has  no  strong 
bias  towards  or  away  from  a sidewall.  The  value  -1.  is 
assigned  if  this  flux  component  (VELO)  is  directed  towards 
the  array  of  detectors.  Most  commonly,  RMCASE  = 0. 

RCSS 

gives  the  depth  of  the  recess,  in  cases  of  a bend  with  a 
trap.  Units  are  the  same  as  for  other  XLTH  values. 

F 

is  an  azimuthal  factor  for  the  radiation  entering  the  bend. 
Normally  it  has  the  value  zero.  (1  - F)  is  applied  as  a 
factor  to  weight  the  component  which  traverses  the  bend 
without  wall  reflections  in  the  bend.  Hence,  for  F=l,  the 
source  entering  the  bend  cannot  exit  without  reflections, 
while  for  F=-l,  the  current  exiting  without  reflections 
(i.e.,  VELO)  is  doubled.  (Note  that  the  configuration  can 
modify  this  component  in  ways  not  expressed  by  basic 
bend  data.) 

RDLNTH 

Measures  the  distance  from  detector  to  the  backscatteri ng 
endwall  of  a straight  duct  segment.  It  is  used  in  RETRN  to 
estimate  and  include  backscattered  current.  It  should  be 
omitted  or  set  zero  if  this  component  is  ignored,  which  is 
the  most  common  case. 

WD 

is  the  width  of  the  entrance  door  into  a room  section 
(R00MX) . 
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WDP 

NW,NL 

AEX 

ARTO 

RAT 

ARAT 


is  the  width  of  the  pattern  of  detectors,  i.e.  the  distance 
from  the  center  of  the  duct  exit  to  the  wall  on  the  side 
where  the  detectors  are  positioned,  in  the  shelter  room 
( ROOMX) . 

(both  fixed  point  numbers  quickly  converted  to  integers)  are 
lateral  and  longitudinal  numbers  of  detectors  in  the  array 
(see  ROOMX).  NW*NL  < 20,  presently,  by  dimension  statement. 

is  a ROOMX  parameter.  It  is  the  area  of  the  entrance  to 
the  room  (IEL=6).  Note  that  IEL  = 7,8  are  alternatives; 
and  DOSX,  rather  than  DOSE  would  be  used  with  ROOMX. 

is  a ratio:  (initial  duct  cross  section  area)/(final  duct 

cross  section  area).  Here,  initial  and  final  refer  to  the 
first  and  last  duct  segments.  For  DOSX,  the  final  aperture 
area  is  the  same  as  AEX. 

is  the  ratio  of  duct  cross-sectional  area  taking  account  of 
a mean  scattering  depth  (PEN)  to  that  with  PEN  = 0.  It 
would  normally  be  unity. 

is  the  ratio  of  cross-sectional  areas  at  a junction  (the 
value  for  the  next  segment  divided  by  that  for  the 
preceding)  if  there  is  a reduction. 
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IV.  Output  Selections  and  Elementary  Examples 


The  input  parameter  I PR  I NT  has  values  0-4  as  al ready  noted.  The 
value  I PR I NT  = 0 suppresses  printout  and  has  no  interest  to  us  here. 

The  value  I PR  I NT  = 4 prints  additional  intermediate  data  and  parameters , 
and  is  useful  mainly  to  inquire  into  the  causes  of  hard-to-expl ai n 
output;  we  will  not  attempt  to  discuss  this  option  further.  Desired 
intermediate  output  may  change  with  one  or  another  particular 
requirement  or  investigation. 

I PR  I NT  = 1 causes  only  input  data,  together  with  minimum  output,  to 
be  printed.  Input  data  have  been  described  in  the  preceding  section  and 
are  readily  recognizable  in  the  printout.  Also  these  cards  and  their 
sequence  can  be  confirmed  from  the  list  INPUT  lines  added  to  the  top  of 
Figs.  2 and  3.  OUTPUT  begins  by  recording  the  input  data  under  suitable 
headings.  This  is  followed  by  values  for  the  dose  and  dose  reduction 
factor  as  indicated  by  printed  titles.  The  heading  OUTPUT  has  been 
added  to  clarify  the  figure  and  does  not  accompany  the  output.  The 
heading  "NEXT  CASE  = .500,"  begins  the  printout  of  results  for  Fig.  2, 
and  this  heading  appears  in  general  with  output  results.  The  number 
.500  gives  distance  down  the  final  straight  segment  specified  by  JVARY, 
and  may  be  zero  if  the  dose  is  computed  at  the  entering  point  for  a new 
segment. 

If  JDEL  = 1,  the  usual  case,  this  distance  will  be  the  seal ed  value 
corresponding  to  XLTH  for  the  segment  identified  by  JVARY.  Thus,  the 
duct  of  Maercker  and  Muckenthaler  (see  Fig.  5)  has  a mean  lateral 
extension  of  /(3‘ x 3 1 ) =3 ' ; and  for  JVARY  = 4 the  value  NEXT  CASE  = .500 
identifies  the  position  .5*3'  = 1.5'  down  the  segment  identified  by  IC  = 
4.  Since  IC  = 2 is  the  first  straight  segment  and  IC  = 3 is  the  first 
bend,  IC  = 4 is  the  second  straight  segment.  The  possibility  for 
incrementing  the  XLTH  value  specified  by  JVARY  by  one  or  more  distances 
DEL  up  to  a total  of  JDEL-1  times  exists.  Then  output  of  the  form 
described  is  repeated,  occurring  JDEL  times  in  succession,  with 
appropriately  changed  headings,  as  shown  in  Fig.  2. 
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INPUT 

2 2031  1.973+000 
15144  1.333+000 
0 


***  M+M,  2-LEG,  THERMAL  FLUX  (USES  BTODR) , 8/29/86  *** 


0 

0.414 

0.000 

0.000 

0.000 

0.000 

2 

3.000 

3.000 

13.500 

0.000 

0.000 

3 

3.000 

3.000 

0.000 

0.000 

0.000 

5 

3.000 

3.000 

1.500 

0.000 

0.000 

7 

0.000 

0.000 

0.000 

0.000 

0.000 

OUTPUT 

•p 

ISORS 

I DOSE  ING 

LPIPE  IEL 

DIMEN 

2 

2 0 

3 1 1. 

973E+00 

IPRBLM 

NEL  IPRINT  JVARY  JDEL  DEL 

1 

5 1 

4 4 1. 

333E+00 

*** 

M,  2-LEG, 

THERMAL  FLUX  (USES 

BTODR),  8/29/86  *** 

IEL 

WIDTH-FT 

HEIGHT-FT 

LENGTH- 

■FT  OTHER 

SELECT 

0 

.414 

.000 

.000 

.000 

.000 

2 

3.000 

3.000 

13.500 

.000 

.000 

3 

3.000 

3.000 

.000 

.000 

.000 

5 

3.000 

3.000 

1.500 

.000 

.000 

7 

.000 

.000 

.000 

.000 

.000 

NEXT  CASE 

.500 

NEXT ( I )=  5 

D0SE=  2. 361E-02 

REDUX  FACTR=  8.459E-03 

NEXT  CASE  1 

.833 

NEXT ( I ) = 5 

D0SE=  4. 369E-03 

REDUX  FACTR=  1.565E-03 

NEXT  CASE  3.166 

NEXT ( I ) = 5 DOSE=  1.215E-03  REDUX  FACTR=  4.355E-04 
NEXT  CASE  4.499 

NEXT ( I ) = 5 DOSE=  5.283E-04  REDUX  FACTR=  1.893E-04 


FIG.  2.  Example  of  elementary  input  and  output  of  DCTDOS,  Data  are  for 
a duct  with  concrete  walls  and  3'x3'  passageway  cross  section. 
Output  is  for  the  2nd  leg,  at  1.5',  5.5',  9.5',  and  13.5' 
beyond  the  bend,  which  in  this  case  is  cubic  in  shape.  The 
(neutron)  source  is  incident  at  45°  at  all  parts  of  the 
entrance  plane  and  all  azimuths. 
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INPUT 


2 

2 0 

3 1 1.973+000 

3 

0 

**3X3X18 

4 1 

2 1 0.000 

1-LEG  DUCT  INTO  ROOM,  THERMAL 

N FLUX, 

8-15-83 

0 

0.000 

0.000  0.000 

0.000 

0.000 

1 

3.000 

3.000  18.000 

0.000 

0.000 

6 

10.000 

8.000  15.000 

3.000 

0.000 

8 

5.000 

1.000  2.000 

9.000 

0.000 

OUTPUT 

$ =$ 

ISORS  IDOSE  ING 

LPIPE  IEL  DIMEN 

2 

2 0 

3 1 1.973E+00 

IPRBLM 

NEL  I PRINT  JVARY  JDEL  DEL 

3 

5 1 

2 1 .00CE+00 

**3X3X18 

1-LEG  DUCT  INTO  ROOM,  THERMAL 

N FLUX, 

8-15-83 

IEL 

WIDTH-FT 

HEIGHT-FT  LENGTH-FT 

OTHER 

SELECT 

0 

.000 

.000  .000 

.000 

.000 

1 

3.000 

3.000  18.000 

.000 

.000 

6 

10.000 

8.000  15.000 

3.000 

.000 

8 

5.000 

1.000  2.500 

9.000 

.000 

7 

.000 

.000  .000 

.000 

.000 

NEXT  CASE  6.000 

DOSE/ENTRANCE  DOSE  W(FT)=  5.00  L(FT)=  15.00 
1 . 0653E-02  9.0436E-03  5.3005E-03 

1.4363E-03  1 . 5541E-03  1 .3522E-03 


DOSE 

4.2036E-02 
5 .6678E-03 


3.5686E-02 
6.1 324E-03 


2.0916E-02 
5 . 3359E-03 


Fig.  3 Example  of  elementary  input  and  output  of  DCTDOS.  A single 
straight  segment  3'x3'xl8'  ends  in  a room  of  dimensions 
lO'xS'xlS1.  The  source  of  thermal  neutrons  is  isotropic  at  all 
parts  of  the  entrance  plane.  See  Fig.  4 for  the  detector 
pattern:  For  6 detectors,  1,  3,  5,  are  opposite  the  entrance 

and  2,  4,  6,  are  against  the  sidewall. 
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DOSX  output  for  I PR  I NT  = 1 gives  arrays  rather  than  single  numbers. 
Thus  in  Fig.  3 output  is  in  the  form  of  two  arrays  --  for  the  dose 
ratio  as  well  as  for  the  dose  itself.  The  heading  for  the  former  is 
DOSE/ENTRANCE  DOSE,  and  for  the  latter,  DOSE,  as  shown  in  Fig.  3.  Both 
arrays  have,  as  their  upper  left  value,  the  entrance  value  to  the  room. 
The  top  line  gives  a sequence  along  the  entering  axis  in  equidistant 
positions  crossing  the  room.  Lower  lines  of  detectors  are  displaced 
laterally  in  the  room  by  equal  amounts  from  entrance  to  sidewall  and 
from  front  to  back,  as  shown  in  Fig.  4 below. 


1 4 7 10 

2 5 8 11 

3 6 9 12 

Fig.  4.  Detector  pattern  for  ROOMX  and  DOSX,  NW  = 3,  NL  = 4. 

This  figure  differs  from  that  for  the  Fig.  3 data  in  its 
proportions  and  in  a non-centered  entrance. 

I PR  I NT  = 2 introduces  several  additional  types  of  information. 

Values  are  printed  for  the  dose  weights;  also  given  are  DOSE  and 
DOSE/ ( ENTRANCE  DOSE)  data  for  the  successive  segment  junctions  of  the 
system.  Thus  to  obtain  the  dose  weights  used,  one  need  run  only  the 
simplest  of  programs  using  I PR  I NT=2 . But  dose  weights  are  also  to  be 
found  in  the  file  PREWTS.  The  types  of  data  which  are  given  with 
IPRINT=2  may  be  required,  or  they  may  provide  a better  understanding  as 
to  why  the  dose  or  dose  reduction  turn  out  to  have  the  values  printed. 

I PR  I NT  = 3 extends  the  output  to  include  spectral  values  together 
with  minimum  intermediate  parameters,  mostly  mean  secants  at  junctions. 

I PR  I NT  - 4 prints  additional  intermediate  values.  Consult  the 
program  listing  in  Appendix  D to  identify  and  possibly  supplement  this 
seldom-used  data  selection. 


16 


V.  Comparison  of  Results  with  Some  Classic  Experiments 


The  comparisons  which  follow  in  this  section  are  absol ute,  based  on 
computations  using  known  source  strengths  and  absolute  detector  data. 

No  relative  comparisons  have  been  included. 

A.  ORNL  Studies.  The  most  extensive  studies  of  neutrons  in  concrete 
ducts  are  those  carried  out  at  ORNL  by  Maerker,  Muckenthal er , and 
others.17-24’43*51  Their  3'  x 3‘  duct  eventually  reached  45'  in  length, 
with  options  of  1,  2,  or  3 bends  as  shown  in  Fig.  5.  In  all  their 
studies  the  source  was  a beam,  directed  through  the  opening  at  an  angle 
of  45°  relative  to  the  duct  axis.  Measurements  were  made  at  points  on 
the  centerline. 

The  ORNL  studies  included  computations  of  the  albedo  angular  and 
energy  functions  for  concrete,  with  a follow-on  computer  Monte  Carlo 
study  based  on  the  resulting  albedo  data.  These  computations  gave  good 
agreement  with  experiments;  and  they  also  permitted  examination  of 
intermediate  spectra  and  orders-of-refl ection  distributions.  Hence 
there  is  a satisfying  completeness  about  the  information  given  in  this 
set  of  reports.  (Note  particularly  Ref.  49.) 

Data  were  developed  for  a thermal  neutron  source  and  for  a fast 
neutron  source.  In  the  latter  case  as  well  as  in  the  former,  thermal 
neutron  flux  was  measured  and  computed.  For  the  latter  case 
measurements  and  computations  for  a Hurst  dosimeter  were  also 
devel oped. 

Fig.  6 shows  the  thermal  neutron  flux  due  to  a thermal  neutron 
source.  Fig.  7 gives  the  thermal  neutron  flux  due  to  an  epithermal 
neutron  source.  Fig.  8 similarly  gives  the  fast  neutron  dose  due  to 
the  "epithermal"  neutron  source.  In  each  case  we  have  added  points 
representing  computations  using  DCTDOS,  and  the  added  points  are 
readily  identified.  The  following  paragraphs  discuss  aspects  of 
these  (DCTDOS)  results. 
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Fig.  5. 


Top  view  Of  source  and  duct 


geometry. 
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THERMAL- NEUTRON  FLUX  (neutrons -cm  / -sec 


ORNL-DWG  66-8671R 


Fig.  6.  Thermal -neutron  flux  in  various  duct  configurations  due  to 
incident  thermal  neutrons 
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SUBCADMIUM  FLUX  (neutrons  cm  * sec' 


0 5 tO  15  20  25  30  35  40  45  50 

CENTER-LINE  DISTANCE  (ft) 


Mg.  7.  Thermal -neutron  fluxes  in  various  duct  configurations  due  to 
incident  epicadmium  neutrons. 
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FAST-NEUTRON  DOSE  RATE  (ergs 


ORNL-DWG  66-1730R3 


£ 


cn 


0 6 12  ,U  18  24  30  36  42 


CENTER-LINE  DISTANCE  FROM  MOUTH  (ft) 


Fig.  8.  Comparison  of  calculated  and  measured  center-line  single- 
collision fast-neutron  dose  rates. 
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A divergence  near  the  source  most  likely  reflects  the  difference 
between  the  experimental  beam  source  and  our  source,  which  covered  the 
aperture  uniformly  and  at  all  azimuths,  and  which  used  interpolated  data 
for  the  45°  obliquity  angle  rather  than  expressing  that  obliquity 
sharply. 

An  increasing  divergence  beyond  22'  in  the  first  straight  segment 
is,  we  believe,  a failure  of  our  model.  But  note  that  a somewhat 
similar  effect  occurs  in  the  ORNL  computations.  It  may  well  indicate 
that  near-grazing  trajectories  are  poorly  described  in  both  models.  In 
our  computations  the  discrepancy  occurs  at  distances  requiring 
extrapolation  beyond  our  data  set;  but  a faulty  extrapolation  formula 
does  not  seem  to  account  fully  for  the  effect. 

A rise  as  a bend  is  approached  is  probably  due  to  back-reflection 
from  surfaces  marking  the  endwall  of  the  bend.  This  can  be  estimated 
if  required. 

High  values  computed  by  DCTOOS  just  following  a bend  can  result  in 
part  from  differences  between  a centerline  detector  and  the  DCTDOS 
detector,  which  averaged  over  the  lateral  plane.  This  will  be  evident 
when  lateral  intensity  variations  are  strong,  as  occurs  when  part 
only  of  the  duct  cross  section  is  exposed  to  intensities  direct  from 
two  segments  preceding  and  possibly  part  of  the  source  plane. 

The  general  agreement  with  the  experiments  appears  to  be 
comparable  with  that  of  the  Monte  Carlo  computations,  if  account  is 
taken  of  the  differences  between  source  and  detector  configurations. 

To  pursue  comparisons  in  greater  detail,  we  computed  both  cumulative 
group  spectra  and  orders-of-refl ection  components  at  several  positions 
in  the  duct.  The  ORNL  studies  published  computed  examples  of  both 
types.  Orders-of-refl ecti on  distributions  were  compared  for  equal 
centerline  distances  both  down  a duct  continuation  and  around  a bend. 
Results  show  that  the  bend  shifts  the  peak  to  values  for  the  orders-of- 
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reflection  which  are  larger  by  4-5,  for  both  first  and  second  bends. 

Our  computations  showed  about  the  same  shift.  But  our  peaks  were 
consistently  located  at  smaller  orders-of-scatteri ng  by  about  the  same 
amount.  We  think  that  this  results  partly  from  our  use  of  a source 
configuration  which  permits  neutron  injection  initially  to  greater 
depths  down  the  first  leg.  But  in  both  ORNL  computations  and  in  ours, 
peak  location  is  not  strongly  sensitive  to  distance  down  a given  leg. 

We  chose  cumulative  rather  than  differential  spectra  for  a 
comparison  because  the  ORNL  energy  group  structure  is  quite  different 
from  ours.49  Fig.  9 gives  the  comparison  between  the  two  types  of 
computed  spectra  at  14'  beyond  the  center  of  the  first  bend.  Below 
about  .05  MeV  the  ORNL  spectrum  rises  higher  to  a maximum  excess  factor 
of  about  4/3  at  an  energy  just  above  the  thermal  component.  This  is 
consistent  with,  although  it  may  not  fully  result  from,  the  larger 
number  of  scatterings  apparently  undergone  by  the  ORNL  flux  at  the  same 
position.  Note  that  the  thermal  components  are  about  equal  at  this 
position,  as  shown  in  Fig.  7 at  29'. 

B.  Armour  Research  Foundation  Experiments  Using  Co60  Gamma  Rays. 

These  experiments  were  not  supplemented  with  computations;  hence 
the  types  of  comparison  are  much  more  limited  than  the  above,  although 
these  data  have  also  be  analyzed  since  by  other  authors. 

In  Fig.  10,  we  show  the  experimental  duct  with  three  legs  used  by 
Terrell  et  al. 70-73  We  have  computed  gamma  attenuation  in  this  concrete 
duct  and  compared  the  data  from  DCTDOS  with  the  experimental  data  of 
Ref.  73.  Our  computations  used  a one-group  approximation  to  the  Co60 
spectrum.  Similarly,  we  do  not  distinguish  between  three  bends  of  the 
Z shape  as  distinct  from  the  U shape,  a difference  which  seems  to  be 
of  minor  importance  here.  Note  that  a poi nt , 50  Curie  source  was 
located  at  S10,  3'  inside  the  first  leg  at  a point  which  we  take  to 
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Cumulative  Fluence  Rate 
Neutrons/cm2/sec/Watt 
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Fig.  9.  Cumulative  epithermal  neutron  fluence  rates. 
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Fig.  10.  Top  view  of  the  6 by  6 foot  concrete  entranceway. 
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be  the  duct  entrance  for  the  DCTDOS  computations.  Recall  that  these 
are  based  on  an  isotropic  input  which  is  assumed  uniform  over  the 
lateral  entrance  plane.  Again  the  centerline  measurements  are  compared 
with  an  average  over  the  lateral  plane  used  by  DCTDOS.  We  have  also 
ignored  source  gammas  back-reflected  into  the  duct. 

The  comparison  is  given  in  Fig.  11.  Our  results  show  less 
satisfactory  agreement  with  the  experiments  than  is  the  case  for 
neutrons,  as  might  be  expected.  The  most  evident  difference  appears 
wthin  and  just  beyond  the  bend,  where  our  computations  do  not 
reproduce  the  feature  that  as  long  as  the  (point-like)  source  can  be 
viewed  by  the  detector,  the  main  fall  off  in  the  second  leg  is 
postponed.  The  consequence  of  this  is  that  our  main  fall  off  occurs  1-2 
feet  later  than  in  the  experiment.  We  believe  here  also  that  at  least 
pert  of  the  explanation  for  this  effect  results  from  a centerline 
approach  vs  an  approach  by  lateral  averaging.  At  the  first  bend,  the 
large  experimental  reduction  occurs  when  the  point  detector  is  first 
shielded  from  the  point  source.  By  contrast,  a plane  source  is  only 
gradually  shielded  from  a plane  detector,  in  this  case  total  shielding 
occurring  nearly  7'  beyond  the  vertical  lines  which  mark  the  turns  to  a 
new  leg  at  the  bend  centers. 

Regarding  the  factor  of  roughly  2 discrepancy  in  the  third  bend,  we 
think  that  this  is  due  to  inadequacies  of  data  expressing  our  model.  In 
addition  to  differences  of  source  and  detector  configuration,  no  attempt 
has  at  this  point  been  made  to  incorporate  and  apply  gamma  ray  albedo 
angular  distributions  to  computations  of  basic  segment  data;  neutron 
albedo  angular  distributions  have  been  used  instead,  as  noted.  The 
substitution  of  more  realistic  data  should  improve  accuracy  somewhat,  at 
the  expense  of  duplicating  much  of  ENTDAT. 


26 


/min 


1000.  0 


100.  0 


10.  0 


1.0  — 


0.  1 


0.  01 
0.  005 

0.  002  - 
0.  001 

0 


A 


& 


Normalized 

here 


O — Point  Source  C 


• - DCTDOS 


6 • 


A 

O 


& 

O 


Centerline  Distance  in  Feet 

_ 1 1 1 1 1 1 L_ 

12  16  20  24  28  32  36 


Z 

shape 


40  44 


Fig.  11.  Measured  gamma  dose  rate  in  6 by  6 foot  concrete  entranceway 
with  two  right  angle  bends  (Z-Shape). 
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VI.  Procedures,  Variables,  and  Basic  Data 


1)  Position,  Energy,  and  Direction  (N,  SO,  SP,  IC,  DLNGTH,  SENVAL). 
We  identify  position  with  the  entrance  or  exit  lateral  planes  in  a 
straight  duct  segment,  except  in  the  special  case  of  a detector  pattern 
in  a room.  _IC_  values  sequence  duct  segments,  and  DLNGTH  is  straight 
segment  length  or  detector  position. 

The  energy  variable  is  expressed  through  the  order  of 
backscatteri ng , _n  or  N=(n+1).  While  there  is  expected  to  be  a 
dependence  of  the  energy  on  the  sequence  of  types  of  backscatteri ng , 
e.g.,  angles  of  incidence  and  emergence,  in  practice  this  dependence  is 
assumed  to  be  weak  and  we  keep  track  only  of  the  number.  In  addition, 
there  is  a dependence  of  the  albedo  on  the  neutron  energy,  which  we 
neglect  because  it  is  not  strong  except  in  special  cases  which  do  not 
concern  us  here.  Our  basic  data  (fcr  concrete)  is  based  on  an  assumed 
albedo  value  of  0.696  throughout  although  more  recent  evaluations 
suggest  that  0.75  is  a oetter  figure.  Rather  than  recompute  the  basic 
data  at  this  point,  we  neglect  the  differences  of  angular  distribution 
as  well  as  spectrum  and  merely  include  a correction  factor 
(0.75/0. 696 ) n to  take  account  of  this,  because  these  differences  do  not 
appear  to  have  a strong  influence  on  segment  outputs.  Since  the  dose 
weights  have  incorporated  the  correct  albedo  factor,  only  the 
correction  factor  ( 1/0 . 696 ) n is  applied  in  DOSE  and  D0SX. 

The  angular  distribution  is  character!' zed , basically,  by  the  ratio 
f 1 ux/current . This  is  an  average  of  the  reciprocal  cosine  over  the 
radiation  current.  (It  is  also  the  reciprocal  of  the  mean 
cosine  averaged  over  the  normalized  flux.)  We  use  this  parameter 
partly  because  it  gives  a simple  way  of  transforming  currents  into 
fluxes  and  vice  versa.  (Such  a relationship  can  only  be  approximate 
because  dose  or  other  energy  weighting  factors  give  variant  flux-to- 
current  ratios. 
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We  have  found  it  convenient  to  use  a derived  quantity 


s 


2 

<sece> 


1 . 


(1) 


Values  of  s for  forward  currents  range  from  -1  to  +1 ; and  the  largest 
value  (<sec9>=l)  corresponds  to  neutrons  travelling  parallel  to  the 
axis  of  penetration,  which  is  nearly  always  the  reference  axis. 
Dependence  on  only  this  single  direction  parameter  is  one  of  the  basic 
limitations  of  our  approach;  but  we  have  found  that  it  works  rather 
well.  It  is  not  difficult  to  identify  configurations  and  directional 
distributions  which  would  be  poorly  described  so  as  to  lead  to 
unsatisfactory  results.  But  such  configurations  are  about  as  likely  to 
suggest  a requirement  for  new  data  types  as  to  suggest  a change  in 
approach. 


In  DCTDOS,  N_  or  NREF  always  identifies  (n  + 1)  rather  than  n_;  and 
S0(_IC_) , SP(JC_)  correspond  to  ^values  for  unreflected  and  reflected 
components,  respectively,  emerging  from  the  segment  J_C.  (But  note  that 
_IC  = 2 is  the  first  real  segment  of  a duct  system  being  evaluated.) 
SENVAL  is  £ value  for  the  current  enteri ng  a bend. 

2)  Dependent  Variables  (VELO,  VELP,  RMTX).  Our  computations  are 
based  on  currents,  rather  than  fluxes,  and  VELO,  VELP  identify  currents 
entering  or  leaving  a duct  segment.  Specifically,  VELP( I C , N ) i s the 
component  of  the  current  of  neutrons  reflected  a total  of  (N-l)  times 
and  emerging  from  the  duct  segment  identified  by  IC,  having  had  at 
least  one  reflection  in  the  IC'th  segment.  Correspondingly,  VEL0( I C , N ) 
is  the  component  of  the  emerging  current  which  passed  through  the  IC'th 
segment  without  any  reflections.  Its  dependence  on  _N  signifies  that  it 
has  no  doubt  been  reflected  in  other  segments,  perhaps  many  times. 

RMTX  is  intended  to  stand  for  Yefl ecti on  matrix^'  The  current 
entering  any  duct  segment  will  be  reflected  0,  1,  2,  ...  n times  in 
that  segment.  These  are  in  addition  to,  and  must  be  combined  with, 
previous  reflections.  The  general  procedure  for  doing  this  is  the 
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evaluation  of  a simple  folding  computation.  This  is  cast  in  the  form 
of  a matrix  operation  in  which  the  incoming  currents,  VELO  and  VELP, 
are  treated  as  vectors  multiplied  by  an  elementary  matrix,  RMTX,  which 
expresses  the  effect  of  the  new  segment.  This  product  then  yields  new 
"vectors"  VELO  and  VELP  constituting  the  output  current.  If  one  has  a 
vector  with  components  VQ,  Vx,  V2,  ...  together  with  a matrix  which  is 
zero  to  the  right  of  the  diagonal,  has  the  value  CQ  on  the  diagonal,  C2 
in  the  positions  just  left  of  the  diagonal,  C2  in  all  the  positions 
just  left  of  the  Cx  elements,  etc.,  then  the  product  is 
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Let  the  subscripts  now  stand  for  numbers  of  reflections,  the  V be 

n 

values  for  the  current  entering  an  element,  and  Cn  for  the  probability 

of  emerging  reflections  in  this  element.  Then  the  matrix  to 

the  right  has  a form  suitable  for  the  corresponding  current  emerging 

from  the  element.  This  assumes  that  the  V are  essentially 

n 

interchangeable  upon  entrance.  Sometimes  we  make  this  assumption;  but 
more  often  we  perform  different  computations  for  VELO  and  VELP , since  the 
entering  angular  differences  seem  often  to  be  important  for  these 
components . 

RMTX  always  is  given  the  general  form  of  the  matrix  shown  above,  but 
if  VELO  and  VELP  are  treated  separately,  RMTX  should  more  accurately  be 
computed  separately  for  each,  a refinement  not  in  DCTDOS  presently. 


3)  Data  Forms  for  Structural  Elements  (BEND,  DUCTO,  DUCT1,  DUCT2) . 

S.  Woolf  has  performed  many  computations  for  duct  segments  of  square 
cross  section  and  different  lengths  or  different  positions  in  a very 
long,  straight  duct.76  Other  data  from  these  computations  explore  more 
general  shapes:  rectangular  rather  than  square  cross  section  shapes. 
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Still  other  computations  are  for  bends  of  square  cross  section  but  with 
lengths  1 x 1 (cubical),  1x2,  etc.;  and  recessed  bends  were  examined. 
Study  of  Woolf's  data  leads  to  the  empirical  result  that  successive 
orders  of  reflection  beyond  n = 1 nearly  always  follow  rather  well  the 
simple  rule 


Cn  = A*an,  for  n > 1 , (3) 

where  Cn  is  the  n'th  reflected  current  at  the  duct  exit,  say.  We 
assign  to  all  of  these  components  combined  the  angular  parameter  P. 
The  n = 0 case  is  usually  quite  different,  and  hence  is  assigned  the 
different  value  SO,  as  already  stated. 

Thus,  to  render  Woolf's  data  in  a form  suitable  for  our 
computations,  we  require  five  parameters:* 

KX  = 1:  Co 

2:  A 

3:  a 

4:  SO 

5:  SP 


*These  types  of  data  are  used  in  the  basic  programs  OCTISO,  DCTSLT,  and 
TURN,  where  the  sequence  of  orders-of-ref 1 ecti on  computed  is  not  that 
of  Eq.  (3)  above,  but  a somewhat  more  complicated  form  which  agrees 
better  with  the  basic  data.  (Note  the  remarks  about  the  subroutine 
ORDST,  in  Appendix  A. 


Fixing  the  square  cross  section  as  basic,  and  the  lxl  bend  shape 
likewise  as  basic,  we  have,  in  the  case  of  straight  ducts,  either  the 
duct  length  or  the  distance  to  the  detector  as  variable:  All  the  above 
types  of  parameter  are  thus  functions  of  this  distance  parameter.  By 
contrast,  in  the  case  of  the  duct  bend  the  dimensions  are  fixed  but 
the  angular  input  is  not;  and  all  of  the  above  types  of  parameter  are 


31 


functions  of  the  value  of  the  incoming  current.  Thus,  straight  duct 
data  are  tabulated  as  a function  of  length  or  distance  to  the  detector, 
both  identified  by  DLNGTH.  Bend  data  are  rendered  as  a function  of  the 
s_  value  (SENVAL)  for  the  entering  current,  which  covers  the  entrance 
surface. 

Straight  duct  data  are  based  on  an  entering  current  proportional  to 
cose  relative  to  the  duct  axis  in  DUCTO,  and  on  entering  currents  with 
mono-obliquity  angles  cose  = 0.25  and  0.75,  in  DUCT1  and  DUCT2.  Other 
data  are  available,  but  we  have  preferred  to  interpolate  or  combine 
these,  as  required.  In  addition,  data  for  finite  ducts  (LP  = 1),  for  a. 
fixed  distance  in  a semi -infinitely  long  duct  (LP^  = 3),  and  for  the 
forward  current  only  in  the  latter  case  (LP  = 2)  are  included  in  ENTOAT. 
LPIPE  and  _LP  specify  the  input  data  type  to  be  used. 

4)  Interpolation,  Extrapolation  of  Segment  Data.  We  rely  mainly  on  a 

simple  2-point  interpolation  procedure  applied  to  the  logarithm  of  data 

points.  Thus,  if  we  wish  data  for  duct  of  length  L which  falls  between 
> * 

tabulated  values  l]  and  L?,  and  if  the  tabulated  values  for  these  two 
points  are  yL  and  y2  respectively,  the  value  y or  the  desired  value  at  L. 
is  given  by 

y = y i (y  2/y  i ) ^ L=”L  1 ^ ^ 1-2-1 1 ^ 

This  is  exact  for  exponential  distributions. 

In  using  such  a simple  interpolation  procedure,  various  data  sets 
are  modified  to  take  account  of  the  following  problems: 

a)  reduce  range  of  variation  by  known  divisor  functions  for  direct 
use  of  DATA  statements, 

b)  ensure  that  neither  yL  nor  y2  is  zero  in  a tabulation, 

c)  counter  a strongly  non-exponential  variation. 

Such  modifications  are  removed  after  an  interpolation. 

Also,  an  extrapolation  procedure  has  been  included  to  make  possible 
rough  estimates  for  extreme  cases  such  as  very  long  and  narrow  ducts. 
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This  is  based  on  the  assumption  that  the  distribution  being  extrapolated 
can  be  described  by  the  general  form 


y = Ax“p  + Bx_m"£  (5) 

for  very  large  values  of  x.  (See  the  discussion  of  the  subroutine 
XTRAP.) 


5)  Combination  of  Distributions.  If  one  wishes  to  combine  additive 
distributions  characteri zed  by  different  values  of  s such  as  VELO  and 
VELP,  a value  for  £ for  the  combination  must  be  obtained  by  taking 
proper  account  of  the  nature  of  s^  as  a mean  value.  Let  us  designate 
two  currents  with  values  s_:  and  s_2  as  ^(s_ L ,cose)  and  C2(s_2,cose) . 
Total  currents  are  designated  c^,  £2,  respectively.  Then 


<sece> 


i 


l 


d(cose)sece[C1(s1  ,cose)  + C2(s2,cose)] 


c + c 
**1  2 


c1*<sec91>  c2*<sec92> 


ci  + c 2 


Cx  + C2 


(6) 


Equation  (1)  can  be  used  to  express  this  equation  in  terms  of  s and  s. , 
i=l ,2. 


6)  Correction  for  Eccentric  Duct  Cross  Sections.  Data  due  to 
Woolf76  show  that  for  the  reflected  component  (VELP)  in  a straight 
duct,  a rough  correction  for  eccentricity  (W  not  equal  to  H)  can  be 
made  which  is  relatively  independent  of  duct  length  over  the  range  of 
our  data.  This  correction  is  not  highly  accurate,  but  it  only  becomes 
large  when  the  ratio  of  side  widths  becomes  more  extreme  than  0.1.  A 
simple  fit  to  this  correction  factor  which  is  independent  of  n has  been 
i ncl uded. 


For  the  unreflected  component  (VELO),  the  correction  is  more 
complex.  It  turns  out  that  duct  length  is  more  important  for  this  case. 
The  empirical  formula  which  we  have  incorporated  attempts  to  modify  the 
simple  form  for  the  unreflected  component  to  include  the  distance 
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parameter  in  about  the  simplest  manner  that  can  represent  the  data. 

The  results  are  not  as  dependable  as  in  the  case  of  the  reflected 
component;  but  the  correction  is  smaller  and  again  it  does  not  become  a 
major  factor  for  duct  eccentricities  much  less  extreme  than  0.1. 

If  e<l  is  the  ratio  of  the  smaller  to  the  larger  of  the  two  lateral 
dimensions,  then  our  correction  factors  have  the  forms 

Kd  = exp( - .109(1 ne) 2)  (7) 

K0  = exp(  - ,06(  1 ne)  2)L/  {V  ) 

The  factor  Kp  multiplies  the  output  list  VELP,  while  K0  multiplies 
the  output  list  VELO,  both  in  DCTISO  and  in  DCTSLT. 

No  correction  of  this  type  is  presently  applied  to  TURN. 

7)  Sudden  Changes  of  Duct  Size.  If  the  duct  cross  sectional  area 
should  suddenly  increase,  the  total  current  would  remain  unchanged  and 
would  rather  quickly  expand  to  fill  the  larger  lateral  space.  By 
contrast,  a sudden  reduction  in  the  cross  sectional  area  of  a duct,  say 
by  a factor  DELR,  would  in  first  approximation  decrease  the  total 
current  by  DELR  but  would  not  otherwise  be  expected  to  change  spectrum 
or  angular  distribution. 

8)  Computations  of  Spectrum  and  Dose.  At  the  conclusion  of  a 
sequence  of  duct  elements  the  currents  VELO  and  VELP  are  represented  by 
angular  parameters  SO,  SP,  respectively  and  by  values  for  all  orders  of 
reflection,  and  the  problem  then  is  to  turn  these  data  into  spectra, 
dose  reduction  factors,  and  possibly  fluences.  Let  us  consider  the 
problem  of  spectra  first. 

A basic  requirement  at  this  point  is  an  evaluation  of  the  neutron 
spectrum  as  a function  of  initial  energy  and  order  of  reflection. 

To  this  end  we  have  defined  a special  energy-loss  problem  in  which 
neutrons  in  a succession  of  23  energy  groups  from  15  MeV  to  thermal  are 
injected  into  the  space  between  plane  concrete  walls  with  directions 
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proportional  initially  to  cose.  If  the  initial  energy  group  is 
identified  with  the  index  IE,  and  the  final  energy  group  by  KE,  the 
resulting  function  A(IE,KE,N)  is  to  be  applied  for  most  of  our 
purposes.*  The  final  spectrum  of  C(KE)  of  a current  V(N)  resulting 
from  an  input  source  with  spectrum  S(IE)  would  then  be  simply 

CO  00 

C(KE)  = l S( IE)  I V(N)*A( IE.KE.N).  (8) 

I E= KE  N=1 

If  the  dose  assigned  to  each  spectral  group  is  designated  D(KE),  it  may 
be  obtained  as  usual  with 


Dose  = l D( KE)  C(KE)  , (9) 

KE=1 


= l S( IE)  l V(N)  l ( D( KE ) A( IE ,KE ,N)  . 

IE=1  N=1  KE=1 


(9’) 


In  eq  (9'),  the  brackets  identify  a two-index  quantity  which  we 
designate  WTSD.  Rearrangement  of  (9')  leads  to  an  alternative  form 
with  an  interchange  between  positions  of  D and  S.  The  corresponding 
bracket  term  is  then  designated  SWTS. 

9)  Flux  and  Flux  Ratios.  The  entering  current  can  readily  be 
weighted  by  dose  factors  to  obtain  a quantity  which  we  refer  to  as 
"current  dose".  Note  that  N = 1 is  the  only  component,  and  presumably 
has  the  the  input  spectrum.  Similarly  one  can  compute  a "current  dose" 
for  other  positions,  as  well  as  the  ratios  to  that  at  the  entrance. 

The  problem  of  the  flux  and  flux  reduction  factors  is  more  complex 
in  that  currents  with  different  angular  distributions  must  be  combined. 


*Later  considerations  suggest  that  it  may  be  desirable  to  supplement 
with  a function  for  mono-obliquity  incidence  cases. 
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Our  approach  is  the  simple  procedure  of  computing  current  dose  for 
emerging  SO  and  SP  components  which  are  then  transformed  by  the  simple 
rul  e 


*current  = <sec0>*current 


(10) 


Appropriate  combinations  and  ratios  then  readily  follow  which  have  an 
accuracy  which  we  think  is  consistent  with  many  applications. 


Additional  information  is  given  in  the  Appendices. 

A listing  of  the  computer  program  is  presented  in  Appendix  D. 
Appendix  A describes  the  various  subroutines  and  their  purpose. 

Appendix  B describes  various  types  of  data  used  and  gives  some 
information  about  data  sources.  Appendix  C gives  additional  examples 
of  input  and  output  of  the  program. 

Albedo  data  has  not  been  printed  here,  partly  because  it  nearly 
doubles  the  number  of  pages  and  partly  because  it  is  most  practical  in 
magnetic  disk  form  for  direct  use  with  computers. 

Despite  its  apparent  promise  and  potential  usefulness,  we  do  not 
regard  the  DCTDOS  program  to  be  completed  in  the  sense  of  some  of 
the  best-known  Monte  Carlo  programs  currently  available.  But  we  hope 
that  the  ideas  incorporated  will  stimulate  additional  work  of  this  type 
leading  to  improvements  and  possibly  applications  to  similar  problems 
in  other  technical  fields  such  as  reactor  engineering. 

In  regard  to  ideas  and  data  developed  herein,  the  author  appreciates 
the  many  contributions  from,  and  interest  expressed  by  Dr.  Stanley 
Woolf,  of  ARCON  Corp.  in  Wakefield,  Massachusetts.  The  use  of  "we"  and 
"our"  throughout  this  manuscript  expresses  the  author's  feeling  of 
partnership  during  the  investigation.  Dr.  G.  L.  Simmons  of  SAIC  in 
San  Diego  assisted  similarly  with  important  aspects  of  the  work. 


VII.  Concluding  Remarks  and  Acknowledgment 
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APPENDIX  A.  Some  Notes  on  Subprograms 


Figure  1 identifies  the  subprograms  of  DCTDOS,  as  well  as  indicating 
the  flow  of  control.  The  list  below  may  also  be  useful  for  quick 
reference.  Following  the  list  are  brief  descriptions. 


a)  Control  Subprogram: 

DCTDOS 

b)  Data  Subprograms: 

~ ENTDAT  (storage) 

DATDCT  (input  data  preparation  and  scaling 
SRSDET  (source  and  detector  data  preparation) 

SND  (preparation  of  source  and  detector  weight  lists  and 
arrays  from  basic  neutron  data) 

SGD  (preparation  of  weight  lists  and  arrays  from  basic  gamma 
ray  data) 

WTLISTS  ( storage  of  computed  weights  in  a permanent  file  called 
PREWTS,  to  avoid  recomputation  with  each  new  run) 

c)  Duct  Segment  Computations: 

DCTISO  (isotropic  source  entering  a straight  duct  segment) 

DCTSLT  (mono-obliquity  source  entering  a straight  duct  segment) 
TURN  (bend  segment  for  different  source  angular  parameters) 
BTODR  (straight  segment  following  a bend) 

ROOMX  (duct  input,  detector  pattern  in  a room) 

d)  Dose  Computations: 

DOSX  (dose  computations  for  ROOMX,  i.e.,  in  a room) 

DOSE  (dose  computations  in,  or  at  the  end  of  a duct) 

e)  Auxiliary  Subprograms: 

RETRN  (makes  estimate  of  backrefl ection  from  end  of  a closed 
duct) 

XTRAP  (a  short  program  for  extrapolation  of  data,  mainly  for 
straight  duct  segments) 

ORDST  (prepares  approximate  orders-of-refl ecti on 
di stributions. ) 

a)  Control  Subprogram  (DCTDOS) . The  pattern  shown  by  figure  1 gives 
a fairly  complete  summary  of  the  options  presented  by  DCTDOS.  For  each 
type  of  case  or  problem  there  is  a loop  indexed  by  IC  which  evaluates  a 
sequence  of  duct  elements.  At  the  conclusion  of  each  problem  or  case 
there  are  options  to  enter  a full  set  of  new  data  (Run  Case),  new 
configuration  data  only  (Problem  Case),  a change  of  a straight  section 
(Length),  or  a new  source  and/or  detector  (Weights). 

In  the  File  #1,  orders  82-85  set  up  a possible  selection  of  reflected 
or  unreflected  component,  although  usually  S=0 . Order  94  applies  LPIPE 
segment  data  to  the  final  segment  only.  Order  118  gives  a new 
application  of  weights  previously  applied  in  the  same  run.  Order  119 
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applies  new  alternative  weights  near  the  end  of  a problem.  Orders  124 
and  148  develop  weights  only  once  for  cases  in  which  a segment  length  is 
varied.  Order  121  develops  gamma  weights  only  once  if  there  are  at  least 
two  bends.  The  orders  from  127  to  132  relate  to  the  computation  of  dose 
for  the  segment  preceding  a room  after  a ROOMX  dose  pattern  computation, 
for  comparisons. 

b)  Data  Subprograms. 

ENTDAT  stores  output  data  for  elementary  duct  segments  (QUCTi  and 
BEND  arrays).  It  also  stores  selected  spectral  source  and  dose 
information  in  most  of  the  remaining  lists.  The  arrays  AWTS  and  AMOD 
represent  obliquity  incidence  data  (sidewall  and  endwal 1 , respecti vely ) 
for  gamma  ray  sources.  These  are  required  by  SRSDET  to  obtain  realistic 
albedo  weights  for  gamma  ray  problems. 

DATDCT  is  mainly  a buffer  program  in  which  input  describing  the 
configuration  is  properly  scaled  and  inserted  in  COMMON  storage.  Input 
is  summarized  in  Table  4,  and  output  of  DATDCT  is  summarized  in  Table 
A5  below.  The  ratio  PEN,  usually  equated  to  zero,  can  be  given  a small 
value  to  take  account  of  the  mean  wall  depth  from  which  reflection 
occurs.  (We  have  made  almost  no  use  of  this  option). 

The  segment  eccentricity  ECC  replaces  the  duct  width  and  height,  and 
COF  is  the  length  scaled  by  the  square  root  of  the  cross  sectional 
area.  DELR  mostly  gives  a cross  sectional  area  ratio  (ARAT)  between 
successive  segments.  But  this  may  not  change,  and  then  0THR=0  so  that 
DELR=1  (order  230).  DELR  is  also  used  to  provide  a required  parameter 
value  for  IEL  = 6 and  8.  DZRO  also  is  given  the  value  of  a required 
parameter  for  IEL  = 3 and  6,  otherwise  parameters  of  occasional 
application  only. 


Table  Al.  Scaled  output  produced  by  DATDCT  and  recorded  with  I PR  I NT  = 2. 

For  IEL  = 6,  ARAT  repeats  the  value  for  the  preceding  segment. 


IEL 

COF 

ECC 

DELR 

DZRO 

SCL 

SO 

PEN 

ZD00R 

RMCASE( I EL= 

=8) 

1 

L//A 

W/H 

ARAT,1 

ARAT 

/RAT 

2 

L//A 

W/H 

ARAT , 1 

ARAT 

/RAT 

3 

L//A 

W/H 

1 

F 

/RAT 

4,5 

L//A 

W/H 

ARAT,1 

ARAT 

/RAT 

6 

W 

H 

L 

WD 

AEX 

7 

RDLNTH//A 

W/H 

ARAT,1 

RDLNTH-.01 

ARTO 

8 

WDP 

NW-1 

NL-1 

ARAT 

ARTO 
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All  the  quantities  named  in  Table  A1  have  the  same  interpretation  as 
in  Table  4,  for  the  same  value  for  IEL.  Note  that  COF  contains  mainly 
scaled  lengths,  ECC  has  eccentricity  ratios,  DELR  has  reduction  factors 
for  constriction. 

SCL  is  a list  (in  COMMON  storage)  of  scale  factors.  Little  use  has 
been  made  of  /RAT,  which  is  ordinarily  unity.  ARTO  is  the  ratio  of 
entrance  area  to  exit  area;  if  there  has  been  a change,  either  single 
or  progressive,  this  scale-factor  can  indicate  the  magnitude  of  the 
total  effect.  But  the  effect  is  computed  segment  by  segment  using  ARAT 
val ues. 

A constriction  confined  to  a segment  entrance  only  is  not  really 
accounted  for  in  DCTDOS,  although  it  would  not  require  much  change. 

Note  that  an  increase  of  segment  cross  section  area  will  spread  the  exit 
flux  over  a larger  area,  while  a reduction  will  reduce  the  total 
entering  flux  without  any  other  effect  on  flux/(exit  area). 

SRSDET  computes  source  and  detector  weights.  It  also  makes  use  of 
three  other  programs:  WTLISTS  provides  for  more  permanent  storage  of  the 
weights,  and  SND  and  SGD  apply  multiple  albedo  data  to  the  computation 
of  orders-of-refl ecti on  data  for  neutrons  and  gamma  rays. 

Orders  274-278  provide  a selection  parameter  for  neutrons,  capture 
gammas,  or  neutrons  plus  capture  gammas. 

Orders  from  280  through  327  are  directed  to  computation  of  orders-of- 
reflection  weights  for  the  gamma-ray  case,  which  depends  strongly  on  the 
configuration  of  the  first  leg  and  first  bend.  The  estimate  here  is 
based  on  use  of  simple  formulae  for  angular  limits  for  a source  at  the 
centerline  point  of  the  entrance  to  a cylindrical  duct  (unreflected  in 
the  first  segment,  the  CUT  limit)  and  the  complementary  CUTP  angular 
range  which  applies  mainly  to  scattered  components  from  the  end-surface 
of  first  leg.  Note  again  that  the  unreflected  from  the  first  segment 
(TVO)  may  be  reflected  following  near-perpendi cul ar  incidence  from  an 
endwall  while  the  reflected  component  (TVP)  exits  the  first  segment 
after  reflection(s)  somewhere  in  the  first  leg,  here  assumed  cylindrical 
for  convenience. 

The  main  task  for  SRSDET  is  the  computation  which  takes  account  of 
spectral  changes  with  multiple  reflections,  and  the  effect  which  this 
has  on  the  multiple-reflection  weights.  This  requires  SND  and  SGD  and 
their  processing  of  the  large  tables  of  multiple-reflection  (albedo) 
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energy  data.  Hence  we  give  below  a more  detailed  analysis  of  SND  and 

SGD,  which  are  the  heart  of  SRSDET. 

SND  Purpose:  to  generate  the  arrays  WTSD,  SWTS,  and  the  list 

ENGAM,  for  neutrons  and  capture  gammas. 

Arguments,  Arrays,  and  Lists: 

SGRUPS:  Input  source  list,  neutrons. 

SWTS:  This  begins  as  the  source  list  which,  multiplied  by  a multi - 

reflection  albedo  energy  array  gives  dose  weight  per  detector 
response  group-strength. 

One  should  note  that  WTSD  and  SWTS  are  compl ementary 
quantities.  After  constructing  each  with  the  orders-of- 
reflection  data,  one  must  still  multiply  the  SWTS  by  the  dose 
weights  in  DGRUPS,  or  correspondingly  the  WTSD  by  the  spectral 
weights  in  SGRUPS,  to  obtain  the  dose  weights  required  (see 
Eqs.9,9'). 

DGRUPS:  Dose  weights  for  the  different  spectral  groups. 

ENGAM:  This  is  a list  of  dose  weights  for  capture  gamma  rays  generated 

in  the  duct  walls.  It  is  assumed  that  the  ratio  of  gamma 
number-current  albedo  to  neutron  number-current  albedo  is 
.095/. 75;  and  it  is  assumed  that  all  the  gammas  are  about  1 
MeV,  with  dose  weight  4.5E-10.  The  indicated  proportion  of 
gammas/neutrons  is  assumed  to  hold  for  each  incident  order  of 
reflection,  with  all  orders  of  reflection  contributing, 
independently  of  spectral  energy.  (These  assumptions  probably 
are  overly  simplistic,  are  untested,  and  should  be  reexamined.) 

ALBEDO:  This  is  the  list  of  emergent  neutron  currents  for  each  order 

of  reflection  (2nd  index)  for  each  energy  group  (first  index), 
for  each  initial  energy  (MEG  value).  These  basic  data  are 
included  into  the  program  by  the  READ  orders. 

Sequence  of  Computations: 

405-416:  Preparatory  values  assigned. 

417-421:  Begin  main  loop,  first  zeroing  ALBEDO,  which  may  not  always 
be  necessary. 

422-426:  Read  the  data  and  assign  the  unref 1 ected  components. 

427-432:  Cumulate  the  contributions  to  both  main  arrays,  plus  ENGAM, 
for  all  source  energies. 
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SGO  Purpose:  for  gammas,  to  generate  the  arrays  WTSD  and  SWTS. 
Arguments,  Arrays,  and  Lists: 

All  are  as  described  in  SND  except  that  ALB  replaces  ALBEDO  of  SND, 
and  that  VEL1  differs  from  ENGAM.  (One  must  interpret  SGRUPS,  DGRUPS, 
and  DWTS,  as  well  as  ALB,  as  applying  to  gamma  rays  in  this  subroutine). 
VEL1:  This  is  a weight  list  for  different  obliquity  angles  with  which 

the  initiating  gamma  rays  may  strike  the  wall.  It  is  used 
because  the  first  order  of  reflection  for  gamma  rays  will  be 
significantly  different,  and  will  affect  all  other  orders  of 
reflection;  hence  some  account  must  be  taken  whether  most  gamma 
rays  first  strike  normally,  grazingly,  or  with  a cosine 
distribution. 

Sequence  of  Computations:  Closely  similar  to  SND. 

451-455:  Zero  the  arrays. 

456-462:  Begin  the  source  energy  loop.  KE  is  the  energy  index  and 
LOB  is  the  obliquity  index.  Note  that  the  energies 
proceed  from  top  down,  each  represented  by  the  five 
obliquity  cases. 

463-467:  Read  the  (mass  storage)  data  and  take  care  of  the 

unreflected  component,  assumed  to  be  unity  and  to  be 
incident  in  proportion  to  VEL1  obliquity  weights. 

468-476:  Cumulate  the  two  main  arrays  over  all  source  energies  and 
obliquities. 

WTLISTS  is  called  by  DCTDOS  at  the  beginning  of  a run  to  set  up  two 
files,  IXFILE  and  WTFILE.  The  latter  contains  precomputed  weights  for 
different  combinations  of  IDOSE,  ISORS,  and  LPIPE,  currently  a maximum 
of  65  sets.  The  former  contains  an  index  which  states  whether  the  set 
is,  in  fact,  stored  for  retrieval  or  should  be  computed.  A second  index 
confirms  the  location  in  the  file,  and  a third  index  gives  the  number  of 
weights.  At  the  end  of  DCTDOS  a final  call  for  WTFILE  stores  the  file 
again  in  PREWTS  including  any  new  weight  lists  that  may  have  been  added. 

This  subroutine  is  called  twice  in  SRSDET.  The  first  call  either 
accesses  precomputed  weights  or  exits  the  routine  because  they  have 
not  as  yet  been  stored.  The  second  call  either  exits  the  routine 
because  they  already  are  included  or  files  a newly  prepared  set  of 
weights  in  PREWTS. 
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c)  Duct  Segment  Computations.  These  constitute  the  main  workiny 
set  for  solutions.  The  three  basic  programs,  DCTISQ,  DCTSLT,  and  TURN 
follow  essentially  the  same  pattern,  namely 

1)  interpolate  required  data  from  tables  in  ENTDAT, 

2)  construct  one  or  more  suitable  backscatteri ng  matrices  for  the 

segment  (see  Eq.  2), 

3)  fold  data  for  the  new  segment  over  the  input  to  obtain  the 

output  as  in  Eq.  2,  and 

4)  modify,  total  (to  VEL)  and  possibly  print  the  results. 

The  program  ROOMX  is  much  more  complex,  because  it  does  not  simply 
focus  on  transmission  of  a current,  but  rather  attempts  to  construct  a 
dose  pattern.  It  treats  separately  the  component  entering  the  room 
which  does  not  reflect  before  contributing  to  the  dose.  (Let  us  call 
this  the  unreflected  component  even  though  it  may  consist  of  many 
earlier  reflections.)  Note  that  we  have  a mean  angular  parameter  for 
each  to  assist  our  description  of  these  two  components. 

Also  given  special  treatment  is  the  once-reflected  component,  which 
is  sensitive  to  the  angular  and  energy  properties  of  the  entering 
current.  All  other  orders  of  reflection  are  lumped  together  and 
treated  with  the  assumption  that  these  contributions  follow  a cosine 
emission  law  from  the  room  surfaces.  The  number  albedo  (.75)  is 
assumed  to  be  modified  (reduced)  by  an  escape  probability  which  is 
taken  to  be  proportional  to  the  ratio  of  entrance  area  to  half  the 
total  wall  area. 

We  believe  that  these  principles  are  reasonably  sound.  Tests 
comparing  results  with  simple  neutron  Monte  Carlo  data  have  been 
promising.  Correspondi ng  tests  for  gamma  ray  data  have  not  been 
carried  out;  and  patterns  for  gamma  ray  cases  should  be  less  accurate. 
Other  ways  to  estimate  these  detector  patterns  can  also  be  devised  and 
could  prove  better. 

BTODR  computes  penetration  down  straight  segments  which  follow 
bends,  and  represents  the  key  to  successful  evaluation  of  mul ti -segment 
ducts.  It  is  a second-generation  program. 
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Consider  a bend  from  the  vantage  point  of  a position  well -down  the 
following  straight  section.  Mainly  visible  is  the  endwal  1 . While 
albedo  from  the  endwal 1 does  not  have  a strictly  isotropic  flux  upon 
emergence,  the  flux  does  not  greatly  differ  from  an  isotropic  distribu- 
tion for  penetration  purposes  down  a straight  section  to  follow.  We 
approximate  this  contribution  with  our  DCTISO  subprogram.  All  three 
other  walls  of  the  bend,  including  the  endwal  1 of  the  first  segment, 
decrease  in  visibility  by  a cosine  factor  with  distance  down  the  next 
segment.  We  approximate  their  contribution  with  our  DCTSLT  program  with 
the  aid  of  our  estimate  of  the  angular  parameters  SO,  SP  of  exit  from 
the  bend  (entrance  into  the  new  segment)  for  components  unscattered  (SO) 
or  scattered  in  the  bend.  It  is  mainly  the  SP  from  the  bend  which  must 
be  apportioned  to  the  walls  which  become  the  new  endwal 1 (for  DCTISO) 
and  the  new  sidewalls.  The  SO  applies  only  to  one  (DCTSLT)  component 
which  happens  to  bypass  the  bend  but  usually  strikes  a wall  near  the  new 
segment  entrance. 

The  remaining  complexities  relate  to  estimates  of  the  number  of 
neutrons  (or  gammas)  which  strike  the  different  walls,  both  reflected 
in  the  bend  and  unreflected,  to  apportion  suitably  their  relative 
contributions  and  mean  angular  incident  at  the  new  segment  entrance. 
Rough  symmetry  between  the  exit  flux  and  that  striking  the  endwal  1 -to-be 
aids  in  making  this  estimate. 

d)  Dose  Computation  Subprograms.  DOSX  is  complicated  by  the  require- 
ment to  evaluate  the  dose  for  the  pattern  used  in  RQOMX.  Otherwise  it 
applies  the  weight  lists  from  SRSDET  to  the  orders  of  reflection 
distributions  at  the  different  specified  positions. 

Both  DOSX  and  DOSE  have  incorporated  several  options,  but  at  present 
the  only  option  used  computes  an  estimated  flux  dose  rather  than  current 
dose  (using  Eq.  10).  Both  also  yield  a reduction  factor  relative  to  the 
initial  position  at  the  duct  entrance. 

e)  Auxiliary  Subprograms.  RETRN  is  of  limited  application,  since  it 
affects  the  final  result  only  in  the  general  vicinity  of  a straight 
segment  endwal 1.  It  is  not  a simple  program,  since  it  must  evaluate 
what  strikes  the  endwall  and  subsequent  back-attenuation.  Possibility 
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of  reduction  due  to  a missing  part  of  the  duct  sidewalls  is  ignored. 

XTRAP  is  a short  subprogram  which  has  taken  different  forms  with 
little  change  in  results  once  the  general  mathematical  form  for 
reduction,  i.e.,  a power  law  down  straight  duct  segments,  is  fixed. 
Present  formula  is  something  of  a hodgepodge  which  works  after  a 
fashi on , 


Y = (AL1  + -^r  ) * 4f  . (Al) 

xu.i  Xxi 

If  the  parameters  of  the  fit  give  AL1<0,  a simpler  form  is  used  instead, 
namely 

Y = ( AL 1 / X X 1 ) (Al1) 

In  the  usual  case,  the  constant  XI  is  pre-assigned  a value  on  the 
basis  of  y a x“XI~*05  for  large  X,  and  AL2,  AL1  are  then  computed. 
Various  methods  for  specifying  the  constant  X_I_  have  been  applied.  This 
present  approach  is  based  on  the  assumption  that  a single  term  is 
sufficient  for  the  large-X  trend.  The  constants  AL1  and  L2^  are  then 
determined  to  fit  the  distribution  well  for  the  largest  tabulated  values 
for  X,  and  give  reasonable  extrapolated  values  at  distances  beyond  the 
tabulation  which  are  not  excessively  large.  The  limiting  case  is  most 
commonly  Y+0  for  X->-«>.  But  a finite  limit  for  Y is  not  ruled  out  and  can 
be  implemented. 

ORDST  stands  for  order-of-refl ection  distribution.  Its  purpose  is 
to  change  a simple  power-law  distribution  fixed  by  two  parameters  into 
a modified  distribution  which  is  more  flexible  and  realistic,  but  which 
can  be  readily  computed  from  the  same  two  parameters.  More  precisely, 
a power  law  can  be  expressed  by  a single  constant  times 

a,  a2,  a3,  a , . .. . , ( A 2 ) 

The  more  flexible  and  applicable  form  which  is  developed  in  this 
subprogram  is  an  overall  factor  times 

a,  a2+ab,  a3+a2b+ab2, (A21) 
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Assignment  of  a value  for  the  parameter  b suggested  by  Woolf  data  for 
elementary  segments76  leads  to  formulae  which  give  _a  and  the  multiplica- 
tive factor  in  terms  of  a.  and  the  original  strength  constant.  Such  a 
change  is  rendered  useful  because  orders-of-refl  ection  distributions  in 
ducts  are  not  simple  power  laws,  and  much  better  descriptions  can  be 
obtained  with  little  more  added  flexibility  than  this  can  give. 


APPENDIX  B.  Sources  of  Data 

a)  S.  Woolf  Duct  Data.  The  unpublished  data  of  ref.  76  constituted  a 
large  number  of  configurations  and  sources,  as  listed  in  Table  B1  below. 
Two  types  of  directional  sources  were  computed,  mono-obliquity  sources 
and  continuous  sources.  The  latter  were  sampled  from  a normalized 
distribution  of  the  form 

f(s,cose)  = 2( 1-s )cos9*d(cose)/[( 1+s) 2-4scose] 3/2  . (Bl) 

The  value  s=0  corresponds  to  the  case  of  incident  isotropic  flux. 

Output  orders-of-refl ection  data  were  given  for  the  zero'th  through 
the  20 ' th  order  of  reflection.  Albedo  data  were  given  similarly.  To 
properly  include  higher  orders  the  Woolf  output  also  gave  values  of  the 
Fourier  Transform, 

oo 

0(p)  = l 0(n)e~np  . ( B2) 

n=I 

These  data  turned  out  to  be  extremely  useful  in  obtaining  simple 
analytic  approximations  to  the  output  distributions-in-n. 

Only  a small  portion  of  the  Woolf  data  were  directly  used  to  obtain 
the  ENTDAT  data  for  DUCT0,  DUCT1,  DUCT2,  and  BEND,  from  which  the  first 
three  parameters  of  Table  4 interpolate.  And  as  noted  in  Section  VI, 
corrections  were  later  necessary  for  an  albedo  of  0.75  instead  of 
0.696. 
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It  was  early  very  clear  that  angular  parameters  would  be  necessary  in 
the  joining  of  different  segments.  But  it  was  not  so  clear  which  or  how 
many  to  use.  To  attempt  perhaps  the  simplest  option  first,  additional 
computations  were  made  to  provide  SO  and  SP,  the  (angular  average) 
parameter  values  listed  in  the  DUCTi  and  BEND  arrays,  for  many  of  the 
most  important  cases. 

The  Table  B1  list  of  cases  were  mostly  recomputed77  in  order  to 
provide  data  for  <sece>,  and  thus  for  SO,  SP. 

Computations  were  also  made  for  eccentric  ducts  with  rectangular 
cross  sections.  The  eccentricity  ratio  e is  used  which  equals 
width/length.  Endpoints  plus  5 internal  planes  of  a duct  of  length  6 
diameters  were  computed.  Eccentricity  ratios  n were  .01,  .25,  .50, 

.75,  and  1.0.  Each  computation  was  performed  for  at  least  6 s values. 
Fortunately  the  sensitivity  to  this  parameter  turned  out  to  be  rather 
low,  except  for  extreme  cases.  This  confirms  also  that  square  cross 
section  data  has  fairly  general  applicability. 


Table  Bl.  Most  of  the  cases  listed  below  were  computed  with  5,000  to 
10,000  histories. 

FINITE,  STRAIGHT,  OPEN  DUCTS,  CONTINUOUS  ANGULAR  INPUT 

8 lengths  (0  to  8 diameters),  each  with  9 input  s values 

FINITE,  STRAIGHT,  OPEN  DUCTS,  M0N0DIRECTI0NAL  ANGULAR  INPUT 
8 lengths  (0  to  6 diameters),  each  with  input  obliquities 
cose.jn=.01,  .25,  .50,  and  .75.  (These  were  redone  to  8 diameters 
with  30,000  histories.) 

SEMI-INFINITE,  STRAIGHT  DUCTS,  CONTINUOUS  ANGULAR  INPUT 

13  depths  (0  to  6 diameters),  each  with  6 s values,  computation 
of  both  flux  and  current. 

OPEN,  STRAIGHT,  FINITE  DUCTS  INTERNALLY,  CONTINUOUS  ANGULAR  INPUT 
Midpoints,  1,  2,  and  4 diameter  lengths,  each  for  6 s values. 

0.5  diameters  depth,  4.5  diameters  length,  each  for  6 s values. 

^ q «.  ii  7 0 11  11  11  11  11  11  11 

CLOSED,  STRAIGHT,  FINITE  DUCTS,  CONTINUOUS  ANGULAR  INPUT 
4,  6 diameters  length,  each  for  6 s values 

OPEN  BEND  SEGMENTS,  SQUARE  CROSS  SECTIONS 

lxl,  2x2,  3x2,  2x1.5,  1.5x2,  7 s values  (continuous  angular  input) 
" " " " " , input  perpendicul ar  to  entrance. 
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b)  SA I Neutron  Multiple-Reflection  Data.38  These  tabulated  data  are 
currently  in  floppy-disk  files.  They  correspond  to  various 
monoenergeti c neutron  fluxes  which  are  backscattered  many  times  between 
two  plane-parallel  concrete  barriers.  At  their  start,  the  angular 
distribution  is  isotropic.  With  each  backscattering , the  emergent 
neutrons  are  classified  into  energy  groups  at  and  lower  than  the  source 
energy,  in  a fairly  standard  set  of  23  energy  groups.  The  concrete 
composition  is  that  referred  to  as  Type  04,  and  has  density  2.34  g/cm3 
and  composition  by  weight  as  shown  in  Table  D3  below: 


Table  B2.  Type  04  concrete  composition  by  fractional  weight. 


ELEMENT 

FRACTIONAL  WEIGHT 

ELEMENT 

FRACTIONAL  WEIGHT 

H 

.0056 

Si 

.3151 

0 

.4979 

S 

.0013 

Na 

.0171 

Ca 

.0829 

Mg 

.0026 

Fe 

.0124 

A1 

.0457 

mul  ti  pi e- 

•albedo  neutron  data 

are  called  and 

used  by  SND,  as 

previously  noted. 


c)  Gamma  Albedo  Data.78  These  data  were  the  basis  for  the  multiple- 
albedo  data  required  for  SGD.  Original  computations  by  Woolf  provided 
detailed  albedo  data  for  single  incidence  only.  A program  was  then 
written  to  apply  these  data  to  the  space  between  plane,  semi -i nf i ni te 
concrete  barriers.  Source  gamma  rays  of  17  different  energy  groups  were 
incident  at  each  of  the  following  obliquity  cosines  (relative  to  the 
perpendi cul ar) , 


cos 0 Q = .05,  .25,  .50,  .75,  1 .0  . 

Albedo  energy  groups  and  emergent  obliquity  cosines  were  then  applied  a 
total  of  10  times  to  evaluate  the  mul ti -reflection  data  in  ALB. 

The  emergent  obliquity  cosines  were  eventually  not  applied  to  the 
gamma  ray  cases,  with  reliance  instead  on  a similarity  between  low 
energy  gamma  and  neutron  albedo  angular  distributions. 
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Our  approach  to  computations  for  the  different  duct  seyment  data 
types  requires  the  application  of  approximate  albedo  anyular  distribu- 
tions in  specialized  Monte  Carlo  studies,  plus  evaluation  of  emeryent 
current-dose  anyular  distributions  characterized  by  the  SO  and  SP 
parameters.  These  permit  ready  estimates  of  fluxes  from  currents  by  the 
simple  rule  previously  given  of 


Flux 


fl  ux 
current 


* current 


a* 


2 

1 + S 


* current  . 


(10) 


Thus  the  emergent  flux  and  flux/current  ratio  are  the  main  angular  data 
recorded,  and  are  used  to  link  the  various  duct  components  together. 


APPENDIX  C.  Remarks  on  Inputs  and  Outputs  with  Examples 

Figure  Cl  gives  both  input  and  output  data  for  a run  with  I PR  I NT =3 . 

To  save  space  let  us  note  that  if  we  had  used  I PR  I NT=2  the  first  page 
would  have  been  identical,  but  would  have  been  followed  by  the  last  five 
lines  of  the  2nd  page  of  these  output  data  to  conclude  the  run. 

The  second  page  gives  source  listing  and  spectrum  at  the  detector,  as 
well  as  a breakdown  of  tne  DOSE  value  into  contributions  from  different 
orders  of  reflection,  termed  DFLXN.  VELO  and  VELP  differ  in  that  the 
former  is  unreflected  in  the  final  segment  while  the  latter  is  reflected 
in  the  final  segment:  both  usually  have  multiple  orders  of  reflection 

in  earl i er  segments.  (Note  that  the  first  column  gives  both  an  energy 
group  list  and  a list  of  orders  of  reflection  (N=n+1). 

The  four  VEL  lists  correspond  to  input  and/or  output  for  successive 
duct  segments.  SO  and  SP  values  below  give  the  angular  distribution 
parameters  _S  for  the  VELO  and  VELP  components,  respectively,  which  are, 
however,  recorded  here  only  for  the  last  VEL  list.  Note  that  for 
forward  currents,  1 < S < -1,  with  the  value  +1  correspond! ng  to  a 
current  focused  along  the  axis,  and  with  the  value  -1  correspond! ng  to 
complete  focus  sideways , i.e.,  90°  from  the  forward  direction.  The 
value  .000  here  in  the  IC=1  column  has  no  meaning,  being  a default  only. 
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INPUT 


32033  1.000+000 
35342  3.000+000 
0 


***  M-M, 

2-LEG, 

EPITH/TH 

ORDERS  OF 

REFL.,  9/2/86 

*** 

0 

0.414 

0.000 

0.000 

0.000 

0.000 

2 

3.000 

3.000 

13.500 

0.000 

0.000 

3 

3.000 

3.000 

0.000 

0.000 

0.000 

5 

3.000 

3.000 

3.500 

0.000 

0.000 

7 

0.000 

0.000 

0.000 

0.000 

0.000 

OUTPUT 

$========$ 

ISORS  IDOSE  ING  LPIPE  IEL  DIMEN 
32033  l.OOOE+OO 
IPRBLM  NEL  IPRINT  JVARY  JDEL  DEL 
35342  3.000E+00 
0 


***  M-M,  2-LEG,  EPITH/TH  ORDERS  OF  REFL.,  9/2/86  *** 


IEL 

WIDTH-FT 

HEIGHT-FT 

LENGTH- 

FT  OTHER 

SELECT 

0 

.414 

.000 

.000 

.000 

.000 

2 

3.000 

3.000 

13.500 

.000 

.000 

3 

3.000 

3.000 

.000 

.000 

.000 

5 

3.000 

3.000 

3.500 

.000 

.000 

7 

.000 

.000 

.000 

.000 

.000 

IEL 

COF 

ECC 

DELR 

DZRO 

SCL 

0 

.41400 

.00000 

.00000 

.00050 

.00000 

2 

4.50000 

1.00000 

1.00000 

.00000 

1.00000 

3 

.00000 

1.00000 

1.00000 

.00000 

1.00000 

5 

1.16667 

1.00000 

1.00000 

.00000 

1.00000 

7 

.00000 

1.00000 

1.00000 

-.01000 

1.00000 

NEXT  CASE 

1.167 

DWTS//SWTD//ENGAM 


4 . 124E+00 

2.965E+00 

2.252E+00 

1.739E+00 

1.320E+00 

9.968E-01 

7.457E-01 

5.576E-01 

4.178E-01 

3 . 131E-01 

2 . 334E-01 

1.747E-01 

1.302E-01 

9 . 715E-02 

7 . 214E-02 

5.317E-02 

3.968E-02 

2.941E-02 

2.207E-02 

1 . 642E-02 

1.214E-02 

4 . 124E+00 

2.965E+00 

2.252E+00 

1.739E+00 

1.320E+00 

9.968E-01 

7.457E-01 

5.576E-01 

4. 178E-01 

3. 131E-01 

2.334E-01 

1.747E-01 

1 . 302E-01 

9.715E-02 

7.214E-02 

5.317E-02 

3.968E-02 

2.941E-02 

2.207E-02 

1.642E-02 

1.214E-02 

.OOOE+OO 

1.690E-10 

1.283E-10 

9.914E-11 

7.523E-11 

5.682E-11 

4.250E-11 

3.178E-11 

2.381E-11 

1.785E-11 

1.331E-11 

9.956E-12 

7.421E-12 

5.537E-12 

4 . 112E-12 

3.031E-12 

2.262E-12 

1.677E-12 

1.258E-12 

9.360E-13 

6.919E-13 

NEXT ( I ) = 
NEXT ( I ) = 
NEXT ( I ) = 
NEXT ( I ) = 


0 DOSE= 

2 DOSE= 

3 DOSE= 
5 DOSE= 


5.833E+00 

8.977E-02 

4.889E-02 

2.347E-02 


REDUX  FACTR=  l.OOOE+OO 
REDUX  FACTR=  1.539E-02 
REDUX  FACTR=  8.381E-03 
REDUX  FACTR=  4.024E-03 


Fig.  Cl.  Example  of  data  from  IPRINT=3. 
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Fig.  Cl 

(Continued) 

ARTO= 

l.OOOE+OO 

MEG.N 

SOURCE 

SPCTRM 

DFLXN 

VELO 

VELP 

1 

€ 

.OOOE+OO 

.OOOE+OO 

.OOOE+OO 

.OOOE+OO 

.OOOE+OO 

2 

1 

.OOOE+OO 

.OOOE+OO 

.OOOE+OO 

.OOOE+OO 

.OOOE+OO 

3 

2. 

, 150E-02 

2.010E-08 

2.365E-04 

5 .648E-05 

5.605E-05 

4 

8, 

.000E-02 

4.532E-07 

4.912E-04 

8.760E-05 

1 .226E-04 

5 

1.405E-01 

1.538E-06 

6.543E-04 

9.016E-05 

1.640E-04 

6 

1. 

, 341E-01 

4.393E-06 

7.069E-04 

7.779E-05 

1 .730E-04 

7 

2. 

. 67  3E-01 

1.074E-05 

6.775E-04 

6 . 113E-05 

1.594E-04 

8 

2. 

. 646E-01 

1.986E-05 

6.023E-04 

4.550E-05 

1.349E-04 

9 

9. 

.440E-02 

4.538E-06 

5.095E-04 

3.277E-05 

1 .079E-04 

10 

3. 

.038E-01 

1.7  57E-05 

4. 166E-04 

2.312E-05 

8.306E-05 

11 

4. 

.430E-01 

6.668E-05 

3.327E-04 

1 .610E-05 

6.227E-05 

12 

3. 

. 922E-01 

2.173E-04 

2.613E-04 

1 . 111E-05 

4.581E-05 

13 

4. 

. 468E-01 

2.341E-04 

2.028E-04 

7.630E-06 

3.325E-05 

14 

2. 

.464E-01 

1.894E-04 

1 . 561E-04 

5.220E-06 

2.389E-05 

15 

2. 

. 350E-01 

4.362E-04 

1.193E-04 

3.563E-06 

1 . 705E-05 

16 

2. 

, 107E-01 

4.200E-04 

9.086E-05 

2.429E-06 

1.211E-05 

17 

2. 

. 099  E -01 

4.961E-04 

6.897E-05 

1.654E-06 

8.569E-06 

18 

1. 

. 492E-01 

3.844E-04 

5.223E-05 

1 . 125E-06 

6.048E-06 

19 

1. 

. 192E-01 

3.246E-04 

3.950E-05 

7.655E-07 

4.261E-06 

20 

1. 

. 497  E -0 1 

4.463E-04 

2.984E-05 

5.206E-07 

2.999E-06 

21 

1. 

.202E-01 

3.684E-04 

1.609E-05 

1.346E-07 

1.586E-06 

22 

9. 

, 570E-02 

3.763E-04 

23 

« 

.OOOE+OO 

1 . 946E-02 

VEL(N)/IC=1 

2 

3 

4 

5 

1 

1. 

, OOOE+OO 

.OOOE+OO 

.OOOE+OO 

.OOOE+OO 

2 

« 

.OOOE+OO 

1 . 927  E -03 

4.540E-04 

.OOOE+OO 

3 

< 

.OOOE+OO 

2. 163E-03 

7.2 1 3E-04 

1 . 125E-04 

4 

* 

.OOOE+OO 

1.849E-03 

7.729E-04 

2.102E-04 

5 

t 

.OOOE+OO 

1.424E-03 

6.869E-04 

2.541E-04 

6 

* 

.OOOE+OO 

1.042E-03 

5 .506E-04 

2.508E-04 

7 

4 

.OOOE+OO 

7.414E-04 

4. 152E-04 

2.205E-04 

8 

« 

.OOOE+OO 

5.185E-04 

3.016E-04 

1.803E-04 

.9 

4 

.OOOE+OO 

3.589E-04 

2. 139E-04 

1.406E-04 

10 

4 

.OOOE+OO 

2.468E-04 

1.495E-04 

1.062E-04 

11 

4 

.OOOE+OO 

1.690E-04 

1.035E-04 

7.837E-05 

12 

4 

.OOOE+OO 

1 . 154E-04 

7 . 114E-05 

5.692E-05 

13 

4 

.OOOE+OO 

7.870E-05 

4.873E-05 

4.088E-05 

14 

4 

.OOOE+OO 

5.360E-05 

3.328E-05 

2.911E-05 

15 

4 

.OOOE+OO 

3.648E-05 

2.269E-05 

2.062E-05 

16 

4 

.OOOE+OO 

2.482E-05 

1.546E-05 

1.454E-05 

17 

.000E+00 

1.688E-05 

1 . 052E-05 

1.022E-05 

18 

4 

.OOOE+OO 

1.148E-05 

7 .157E-06 

7.173E-06 

19 

4 

.OOOE+OO 

7.802E-06 

4.868E-06 

5.027E-06 

20 

4 

.OOOE+OO 

5.304E-06 

3.310E-06 

3.519E-06 

21 

4 

.OOOE+OO 

3.606E-06 

2.250E-06 

1 . 721E-06 

S0= 

4. 

. 140E-01 

4 . 547  E -01 

2.680E-01 

3.602E-01 

SP= 

4 

.OOOE+OO 

3. 532E-01 

.OOOE+OO- 

-7.022E-02 

NEXT  CASE 

4.167 

NEXT) 

;i; 

1=  0 DOSE=  5.833E+00  REDUX 

FACTR=  l.( 

100E+00 

NEXTl 

i 

1=  2 DOSE=  8.977E- 

-02  REDUX 

FACTR=  1.! 

539E-02 

NEXT! 

i 

1=  3 DOSE=  4.889E- 

-02  REDUX 

FACTR=  3.: 

381 E -03 

NEXTl 

[i; 

1=  5 OOSE=  1.410E- 

-03  REDUX 

FACTR=  2.418E-04 
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Fig.  Cl  (Continued) 

ART0= 

= 1.000E+00 

MEG,N  SOURCE 

SPCTRM 

DFLXN 

VELO 

VELP 

1 

.000E+00 

.OOOE+OO 

.OOOE+OO 

.OOOE+OO 

.OOOE+OO 

2 

.OOOE+OO 

.OOOE+OO 

.OOOE+OO 

.OOOE+OO 

.OOOE+OO 

3 

2.150E-02 

1.234E-09 

1.472E-05 

9 .337E-06 

1 .242E-06 

4 

8.000E-02 

2.641E-08 

2. 598E-05 

1.277E-05 

3.863E-06 

5 

1.405E-01 

9.033E-08 

3.287E-05 

1.236E-05 

6 .445E-06 

6 

1.341E-01 

2.550E-07 

3.557E-05 

1.031E-05 

7.994E-06 

7 

2.673E-01 

6.177E-07 

3.51  IE-05 

7.945E-06 

8.410E-06 

8 

2.646E-01 

1 . 131E-06 

3.268E-05 

5.842E-06 

8.0Q3E-06 

9 

9.440E-02 

2.576E-07 

2.924E-05 

4 . 176E-06 

7.136E-06 

10 

3.038E-01 

9.958E-07 

2.547E-05 

2.931E-06 

6.088E-06 

11 

4.430E-01 

3.728E-06 

2.175E-05 

2.035E-06 

5.036E-06 

12 

3.922E-01 

1 . 187E-05 

1.832E-05 

1.402E-06 

4.073E-06 

13 

4.468E-01 

1.267E-05 

1.528E-05 

9.613E-07 

3.241E-06 

14 

2.464E-01 

1.022E-05 

1.264E-05 

6.571E-07 

2.546E-06 

15 

2.350E-01 

2.334E-05 

1.039E-05 

4.483E-07 

1.981E-06 

16 

2. 107E-01 

2.24QE-05 

8.498E-06 

3.Q55E-07 

1.529E-06 

17 

2.099E-01 

2.645E-05 

6.924E-06 

2.079E-07 

1.173E-06 

18 

1.492E-01 

2.049E-05 

5.624E-06 

1.415E-07 

8.950E-07 

19 

1.192E-01 

1.731E-05 

4.555E-06 

9 .623E-08 

6.801E-07 

20 

1.497E-01 

2.378E-05 

3.681E-06 

6.544E-08 

5. 149E-07 

21 

1.2Q2E-01 

1.967E-05 

2.116E-06 

3 .207E-08 

2.768E-07 

22 

9.570E-02 

2.012E-05 

23 

•OOOE+OO 

1.195E-03 

VEL(N)/IC=1 

2 

3 

4 

5 

1 

1.000E+00 

.OOOE+OO 

.OOOE+OO 

.OOOE+OO 

2 

.000E+00 

1 . 927E-03 

4 .540E-04 

.OOOE+OO 

3 

.000E+00 

2.163E-03 

7.213E-04 

1.058E-05 

4 

.000E+00 

1.849E-03 

7.729E-04 

1.664E-05 

5 

.000E+00 

1.424E-03 

6.869E-04 

1 . 88  IE -05 

6 

.000E+00 

1.042E-03 

5.506E-04 

1.831E-05 

7 

.OOOE+OO 

7.414E-04 

4. 152E-04 

1.636E-05 

8 

.000E+00 

5 . 185E-04 

3 .016E-04 

1.385E-05 

9 

.000E+00 

3.589E-04 

2. 139E-04 

1 . 131E-05 

10 

.OOOE+OO 

2.468E-04 

1.495E-04 

9.019E-06 

11 

.OOOE+OO 

1.690E-04 

1.035E-04 

7 .070E-06 

12 

.OOOE+OO 

1 . 154E-04 

7 . 1 14E-05 

5.475E-06 

13 

.OOOE+OO 

7.870E-05 

4.873E-05 

4. 202E-06 

14 

.OOOE+OO 

5.360E-05 

3.328E-05 

3.203E-06 

15 

.OOOE+OO 

3.648E-05 

2.269E-05 

2.429E-06 

16 

.OOOE+OO 

2.482E-05 

1.546E-05 

1.834E-06 

17 

.OOOE+OO 

1.688E-05 

1 . 052E-05 

1.381E-06 

18 

.OOOE+OO 

1 . 148E-05 

7.157E-06 

1.037E-06 

19 

.OOOE+OO 

7.802E-06 

4.868E-06 

7 . 763E-07 

20 

.OOOE+OO 

5 .304E-06 

3 . 31 OE -06 

5.804E-07 

21 

.OOOE+OO 

3.606E-06 

2.250E-06 

3.088E-07 

S0= 

4 . 140E-01 

4.547E-01 

2.680E-01 

7.438E-01 

SP= 

.OOOE+OO 

3.532E-01 

.OOOE+OO 

2.636E-01 
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The  second  case  (4.167)  follows  the  same  printout  pattern. 

Figure  C2  reproduces  the  input  for  all  the  ORNL  neutron  duct 
penetration  sources  of  Figs.  6-9.  Only  a single  run  was  requi red 
because  different  sources  and  detectors  were  cycled  at  each  detector 
point.  On  an  IBM  XT,  with  the  8087  chip  and  storage  of  the  ALBEDO, 

ALB,  PREWTS  files  in  the  RAM  memory,  this  run  required  about  5 minutes. 
It  should  be  noted  that  the  (unit)  thermal  neutron  source  required 
additional  normalizing  of  the  data  by  a factor  of  very  nearly  2 
(neutrons/cm2)  to  obtain  the  experimental  source  strength.  With  use  of 
IPRINT  = 15  the  many  cases  require  about  9 pages  of  output. 

Figure  C3  shows  the  input  for  the  ARF  C060  gamma  ray  experiments. 
Only  a single  detector  response  has  been  included.  Note  that  in  this 
case  the  first  input  data  card  suffices  for  1-leg,  2-leg,  and  3-leg 
cases.  Time  required  was  a little  over  4 minutes.  In  every  computa- 
tion about  5/3  minutes  is  required  at  the  start  to  set  up  the  PREWTS 
data  file  into  the  arrays  consulted  in  SRSDET.  Actual  computation 
time  of  problems  can  be  better  estimated  by  subtracting  this  setup 
interval . 

The  references  include  many  other  interesting  cases,  including 
design  data  such  as  are  in  Refs.  3 and  4,  additional  experimental  data 
from  many  of  the  references  listed,  Monte  Carlo  studies,  and  the 
various  semi -analytic  and  semi -empi ri cal  studies.  Note  that  concrete- 
walled  cavities  are  presently  required,  or  in  the  case  of  gamma  rays 
moderately  low-Z  walls.  The  program  is  also  written  only  for  near- 
right-angle  bends.  Because  of  the  similarity  of  the  137Cs  ARF  results 
to  the  60Co  results  we  have  not  included  Cesium  source  computations 
in  this  report. 

Preliminary  applications  have  also  been  made  to  some  duct  configura- 
tions for  possible  Civil  Defense  structures.  These  have  not  been 
included  in  this  report  because  comparison  data  is  lacking. 
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INPUT 

2 

2 0 

3 1 

1.000+000 

2 

•5 

3 I 

2 11 

1.333+000 

2 

3 4 

0 3 

1.000+000 

2 

3 3 

0 3 

3.600+005 

2 

3 3 

1 3 

3.600+005 

***MAERCKER-MUCKEN,  1-LEG, 

THERMAL  1 

MEUTS,  8/31/86* 

** 

0 

0.414 

0 e 000 

0.000 

0.000 

0.000 

2 

3.000 

3.000 

0.000 

0.000 

0.000 

7 

0.000 

0.000 

0.000 

0.000 

0.000 

2 

2 0 

3 1 

1.000+000 

2 

o 

5 1 

4 8 

1.333+000 

•3 

2 

3 4 

0 3 

1.000+000 

2 

3 3 

0 3 

3.600+005 

2 

3 3 

1 3 

3.600+005 

***M+M, 

2-LEG, 

THERMAL  FLUX,  8-15-1 

33  *** 

0 

0.414 

0.000 

0.000 

0.000 

0.000 

2 

3.000 

3.000 

13.500 

0.000 

0.000 

3 

3.000 

3.000 

0.000 

0.000 

0.000 

5 

3.000 

3.000 

0.500 

0.000 

0.000 

7 

0.000 

0.000 

0.000 

0.000 

0.000 

2 

2 0 

3 1 

1.000+000 

3 

r> 

7 1 

6 6 

1.333+000 

>} 

3 

3 4 

0 3 

1.000+000 

3 

3 3 

0 3 

3.600+005 

3 

3 3 

1 3 

3.600+005 

***MAERCKER-MUCKEN,  3-LEG 

i,  THERMAL 

NEUTS,  7/26/83 

[ *** 

0 

0.414 

0.000 

0.000 

0.000 

0.000 

2 

3.000 

3.000 

13.500 

0.000 

0.000 

3 

3.000 

3.000 

0.000 

0.000 

0.000 

5 

3.000 

3.000 

12.000 

0.000 

0.000 

3 

3.000 

3.000 

0.000 

0.000 

0.000 

5 

3.000 

3.000 

0.500 

0.000 

0.000 

7 

0.000 

0.000 

0.000 

0.000 

0.000 

Fig*  C2*  Input  Data  for  computations  of  the  ORNL  duct. 
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INPUT 

55031  3.319+009 
1 3 1 2 12  3.000-001 


0 


***  ARF 

1-LEG, 

50 

CURIE  C0-60 , RAD 

DOSE,  9/1/86 

*** 

0 

0.000 

0.000 

0.000 

0.000 

0.000 

1 

6.000 

6,000 

1.800 

0.000 

0.000 

7 

-1 

0 

***  ARF 

0.000 
5 1 

0.000  0.000 
4 10  3.000-001 

0.000 

0.000 

2-LEG, 

50 

CURIE  C0-60 , RAD 

DOSE,  9/1/86 

*** 

0 

0.000 

0.000 

0.000 

0.000 

0.000 

1 

6.000 

6.000 

10.000 

0.000 

0.000 

3 

6.000 

6.000 

0.000 

0.000 

0.000 

5 

6.000 

6.000 

1.000 

0.000 

0.000 

7 

-3 

0 

***  ARF 

0.000 
7 1 

0.000  0.000 
6 6 3.000-001 

0.000 

0.000 

3-LEG, 

50 

CURIE  C0-60 , RAD 

DOSE,  9/1/86 

*** 

0 

0.000 

0.000 

0.000 

0.000 

0.000 

1 

6.000 

6.000 

10.000 

0.000 

0.000 

3 

6.000 

6.000 

0.000 

0.000 

0.000 

5 

6.000 

6.000 

8.000 

0.000 

0.000 

3 

6.000 

6.000 

0.000 

0.000 

0.000 

5 

6.000 

6.000 

1.000 

0.000 

0.000 

7 

0.000 

0.000 

0.000 

0.000 

0.000 

Fig.  C3.  Input  Data  for  computations  of  the  ARF  duct. 
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APPENDIX  D.  FORTRAN  Listing  of  DCTDQS 


The  listing  of  DCTDOS  subprograms  which  follows  can  be  assembled 
and  run  using  Microsoft  FORTRAN  software  with  a microcomputer . It 
also  requires  three  additional  files,  called  NBS-TD-G,  NBS-TD-N,  and 
PREWTS.  The  first  of  these  file  contains  gamma  ray  multi-reflection 
albedo  data,  the  second  contains  neutron  mul ti -refl ection  data,  and  the 
third  must  have  a format  acceptable  to  the  subroutine  WTLISTS,  but  is 
initially  full  of  zeros.  All  three  files  are  between  45,000  and  56,000 
bytes  in  size.  Running  times  earlier  quoted  did  not  include  preparatory 
storage  of  the  files  in  RAM  memory. 

For  assembly  purposes  we  have  divided  DCTDOS  '’nto  four  parts  roughly 
equal  in  size  because  of  a 64,000  byte  memory  limit.  After  assembly  of 
each  part  the  object  programs  have  been  linked  into  a run  program  of 
size  roughly  180,000  bytes. 

The  subprograms  have  been  little  affected  by  the  change  from  a 
mainframe  to  a microcomputer  with  the  exception  of  ENTDAT,  which 
required  a special  format. 
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1 : i 

2: 

3: 

4: 

5: 

6: 

7: 

8: 

9: 

10: 

11: 

12: 

13: 

14: 

15: 

16: 

17: 

18: 

19: 

20: 

21: 

22: 

23: 

24: 

25: 

26: 

27: 

28: 

29: 

30: 

31: 

32: 

33: 

34: 

35: 

36: 

37: 

38: 

39: 

40: 

41: 

42: 

43: 

44: 

45: 

46: 

47: 

48: 

49: 

50: 

51: 

52: 

53: 

54: 

55: 

56: 

57: 

58: 

59: 

60: 

61: 

62: 

63: 

64: 

65: 

66: 

67: 

68: 

69: 

70: 

71: 

72: 

73: 

74: 

75: 

76: 

77: 


FILE  1 


DCTDOS  IS  TO  CALCULATE  DOSE  IN  A DUCT-BEND-ROOM  CONFIGURATION 
COMMON  DIAMS ( 18 ) , DUCTO ( 3 . 17 . 6 ) . DUCT1 ( 3 . 17 , 6 ) , DUCT2 (3,17,6) .IEL.LP. 
X BEND (10.6). SGRUPS ( 23 ) , DWTS ( 21 ) . CRRNTO (9.10). RCRRNT (9,10) 

COMMON  RVCT ( 21 ) , RMTX ( 21 , 21 ) . SO ( 9 ) . SP ( 9 ) . DZRO ( 9 ) , NEXT ( 9 ) . COF ( 9 ) . 

X ECC ( 9 ) , DELR ( 9 ) ,SNYNEU(23) .SINVAL(IO) , COl , C02 , VEL (9,21) , IPRINT 
COMMON  RCOSO . RCOSP , RCOS . DGRUPS ( 23 ) , VELO (9,21). VELP (9,21) 

X , VEL1 ( 5 ) .SCALE. SWTS (23. 21) ,WTSD(23.21) 

COMMON  IDOSE , HNDRSN ( 23 ) , SEPI ( 23 ) , STH ( 23 ) , SINR ( 23 ) . THMLN ( 23 ) 

X, ENGAM ( 21 ) , DIMEN, SCOBLT( 23) ,AWTS(6,5) ,SGAM(23) ,DGAM(23) ,AMOD(5,5) 
DIMENSION  SVELO (22.20), SVELP (22.20). SCL ( 9 ) , SELECT ( 9 ) 

DIMENSION  IPR ( 6 ) . ISO ( 6 ) , IDO ( 6 ) . NGI ( 6 ) . ILP ( 6 ) . DIMI ( 6 ) 

1 FORMAT ( // ///20X, 10H$========$) 

66  FORMAT (13,13,13,13,13, 1PE10 . 3 , 30H  ( IPR, ISO . IDO , ING . ILP . DIMEN ) ) 

2 FORMAT (F10. 3) 

3 FORMAT ( F10 . 3 , 7D13 . 3 ) 

4 FORMAT ( 7 ( 1PE10 . 3 ) ) 

5 FORMAT (1H  ) 

6 FORMAT (15, 5F10 . 5) 

7 FORMAT (1H1) 

8 FORMAT ( 2 5H  NEXT  CASE.F7.3) 

9 FORMAT ( 515 , 1PE10 . 3 ) 

10  FORMAT ( 55H  IPRBLM  NEL  IPRINT  JVARY  JDEL  DEL  ) 

11  FORMAT (55H  ISORS  IDOSE  ING  LPIPE  I EL  DIMEN  ) 

12  FORMAT (55H  I EL  COF  ECC  DELR  DZRO  SCL  ) 

OPEN ( 10 , FILE= ' B : NBS-DT-G ' ) 

OPEN ( 11 , FILE= ’ B : NBS-DT-N ' ) 

OPEN ( 05 , FILE= ' ') 

OPEN ( 06 , FILE= ' A : PAPER ' ) 

OPEN ( 07 , FILE= ' B : PREWTS ' ) 

IOX=-l 

CALL  WTLISTS ( JSTT . JXWT , IOX . NMAX . DWTS , SWTD . ENGAM , IDOSE . ISORS . LP ) 

CALL  ENTDAT 

MODG=0 

M0DN=0 

13  READ (5, 9)  ISORS , IDOSE . ING . LPIPE . IEL . DIMEN 
IPRINT=IEL 

IF ( IPRINT. GT.O)  WRITE (6. 7) 

16  IDEL=0 

READ (5.9)  IPRBLM , NEL . IPRINT . JVARY , JDEL , DEL 
IF ( IPRINT. EQ.O)  GO  TO  50 
WRITE (6,1) 

WRITE (6. 11) 

WRITE (6.9)  ISORS , IDOSE , ING . LPIPE , IEL . DIMEN 
WRITE (6. 10) 

WRITE (6.9) IPRBLM . NEL . IPRINT . JVARY . JDEL , DEL 

50  READ (5. 9)  JDOS 
WRITE ( 6 , 9 ) JDOS 
IDOS=0 

IF ( JDOS . LT . 1 ) GO  TO  68 
DO  67  1=1. JDOS 

READ (5,9)  IPR ( I ) , ISO ( I ) , IDO ( I ) . NGI ( I ) , ILP ( I ) , DIMI ( I ) 

IF (IPRINT. EQ.O) GO  TO  69 

WRITE (6.9)  IPR ( I ) , ISO ( I ) , IDO ( I ) , NGI ( I ) , ILP ( I ) . DIMI ( I ) 

69  IPR ( JDOS+1 ) = IPRBLM 
ISO( JDOS+1) =ISORS 
IDO ( JDOS+1 )=IDOSE 
NGI ( JDOS+1 )=ING 
ILP ( JDOS+1 )=LPIPE 
DIMI( JDOS+1) =DIMEN 

67  CONTINUE 

68  CALL  DATDCT(NEL. SCL. SELECT) 

IF (IPRINT.LT. 2)  GO  TO  14 
WRITE (6,12) 

WRITE (6,6)  (NEXT ( I ) ,COF(I) .ECC(I) , DELR(I) ,DZRO(I) .SCL(I) ,1=1, NEL) 

14  IC=1 

IF (IPRINT. EQ.O)  GO  TO  51 
WRITE (6,8)  COF ( JVARY ) 

51  CONTINUE 
VELO (IC, 1) =1 . 

VEL ( IC , 1 ) =1 . 

VELP ( IC , 1 ) =0 . 

RVCT ( 1 ) =1 . 

SO(IC) =COF(l) 

DO  15  1=2.21 
VELO ( IC , I ) =0 . 

VELP(IC.I) =0. 
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78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 


VEL ( IC , I ) =0 . 

15  RVCT ( I ) =0 . 

30  IF (NEL-IC)  1000.1000.31 

31  IC=IC+1 

S= SELECT (1C) 

IF(S.EQ.O.O)  GO  TO  33 
V0=1.+0.5*S*(1.-S) 

VP=1.-0.5*S*(1.+S) 

DO  40  N=1 , 21 

VELO ( IC-1 , N) =VO*VELO (IC-l.N) 

40  VELP ( IC-1. N)=VP*VELP( IC-l.N) 

33  IEL=NEXT ( IC ) 

SCALE=SCL ( IC) 

LP  = 2 

RDLNTH  = DZRO ( NEL ) 

IF ( IC+l-NEL) 35,32,32 

32  LP=LPIPE 
C 

35  GO  TO  (100.110,300,500,500,400,200,210) , IEL 
C 

100  DLNGTH=COF ( IC) 

CALL  DCTISO ( DLNGTH , IC) 

GO  TO  30 
C 

110  DLNGTH=COF ( IC ) 

CALL  DCTSLT ( DLNGTH, IC) 

GO  TO  30 
C 

300  CALL  TURN ( SEN, IC) 

GO  TO  30 
C 

400  CALL  ROOMX ( IC , SVELO , SVELP ) 

GO  TO  30 
C 

500  DLNGTH=COF ( IC ) 

CALL  BTODR ( DLNGTH. IC) 

GO  TO  30 
C 

200  IEL=NEL 

IF (NEXT ( IC-1 ) .EQ.8)  GO  TO  204 
IF ( IPRBLM . LT , 0 ) GO  TO  205 
IF( JDOS . GE. 1)  GO  TO  201 
IF ( IDOSE . NE . 5 ) GO  TO  202 
IF (NEL. GE. 5)  GO  TO  202 

201  CALL  SRSDET ( ISORS , ING , MODG , MODN ) 

GO  TO  205 

202  IF(IDEL.LT.l)  CALL  SRSDET ( ISORS , ING , MODG . MODN ) 
GO  TO  205 

204  IF  ( IDOS  . GT  . 0 ) GO  TO  205 
IC=IC-2 

COF ( IC) =COF ( IC+2 ) 

SCL(IC) =SCL(IC+2) 

DELR(IC) =DELR ( IC+2 ) 

ECC(IC) =ECC ( IC+2 ) 

DZRO ( IC ) =DZRO ( IC+2 ) 

205  CALL  DOSE ( IC ) 

IF( JDOS.LT.l)  GO  TO  220 
IDOS=IDOS+l 
IPRBLM=IPR ( IDOS ) 

ISORS=ISO ( IDOS ) 

IDOSE=IDO ( IDOS ) 

ING=NGI ( IDOS ) 

LP=ILP ( IDOS ) 

DIMEN=DIMI ( IDOS ) 

206  IF ( IDOS. GT. JDOS) GO  TO  220 
IF ( IPRINT . EQ . 0 ) GO  TO  207 

WRITE (6,66)  IPRBLM . ISORS . IDOSE . ING , LP , DIMEN 

207  IF (NEXT (IC-1) .EQ.8)  GO  TO  205 
GO  TO  200 

210  IEL=NEL 

IF(IDEL.LT.l)  CALL  SRSDET ( ISORS . ING , MODG . MODN) 
CALL  DOSX(IC, SVELO, SVELP) 

IF (JDOS.LT.l)  GO  TO  212 
IDOS=IDOS+l 
IPRBLM=IPR ( IDOS ) 

ISORS=ISO ( IDOS ) 

IDOSE = IDO ( IDOS ) 

ING=NGI ( IDOS ) 

LP=ILP ( IDOS ) 

DIMEN=DIMI (IDOS) 

IF ( IPRINT. EQ.O) GO  TO  211 
IF  ( IDOS . GT . JDOS ) GO  TO  212 

WRITE (6.66)  IPRBLM . ISORS . IDOSE . ING . LP . DIMEN 


161 

162 

163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 

181 

182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 

201 

202 

203 

204 

205 

206 

207 

208 

209 

210 

211 

212 

213 

214 

215 

216 

217 

218 

219 

220 

221 

222 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 

241 

242 

243 


FILE  1 


Continued 


211  IF ( IDOS . LE . JDOS ) GO  TO  210 

212  IF(IC.EQ.NEL)  GO  TO  240 
IF(IC.LT.NEL)  IC=IC+1 
IF(IC.EQ.NEL)  GO  TO  233 

220  IDEL=IDEL+1 

IF (IDEL-JDEL)  230.240.240 
230  COF ( JVARY ) = COF ( JVARY ) + DEL 
IF ( JDOS . LT . 1 ) GO  TO  14 
233  IDOS=0 

IF (NEXT ( IC-1 ) .EQ.8)GO  TO  33 
GO  TO  14 

240  IF ( IPRBLM**2-4 ) 16.13.1000 


1000  IOX=l 

CALL  WTLISTS ( JSTT . JSWT , IOX . NMAX . DWTS . SWTD . ENGAM . I DOSE , ISORS . LP ) 

STOP 

END 

SUBROUTINE  DATDCT ( NEL . SCL , SELECT ) 

COMMON  DIAMS ( 18 ) , DUCTO ( 3 . 17 , 6 ) , DUCT1 ( 3 , 17 . 6 ) , DUCT2 (3,17,6) .IEL.LP, 
X BEND (10.6). SGRUPS ( 23 ) . DWTS ( 21 ) . CRRNTO (9,10). RCRRNT (9.10) 

COMMON  RVCT ( 21 ) . RMTX ( 21 , 21 ) . SO ( 9 ) . SP ( 9 ) . DZRO ( 9 ) . NEXT ( 9 ) . COF ( 9 ) . 

X ECC ( 9 ) , DELR ( 9 ) ,SNYNEU(23) .SINVAL(IO) , COl , C02 , VEL (9,21) . IPRINT 
COMMON  RCOSO . RCOSP , RCOS . DGRUPS ( 23 ) , VELO (9,21), VELP (9.21) 

X , VEL1 ( 5 ) .SCALE, SWTS (23, 21) ,WTSD(23,21) 

COMMON  I DOSE , HNDRSN ( 23 ) . SEPI ( 23 ) , STH ( 23 ) . SINR ( 23 ) . THMLN ( 23 ) 

X , ENGAM ( 21 ) , DIMEN , SCOBLT ( 23 ) . AWTS (6.5), SGAM ( 23 ) . DGAM ( 23 ) , AMOD (5.5) 
DIMENSION  SCL ( 9 ) . HGHT ( 9 ) , WDTH ( 9 ) . XLTH ( 9 ) . AREA ( 9 ) , OTHR ( 9 ) , SELECT ( 9 ) 
CHARACTER* 80  TITLE 

2 FORMAT (F10. 3) 

3 FORMAT ( F10 . 3 , 7D13 . 3 ) 

4 FORMAT ( 7 ( 1PE10 . 3 ) ) 

5 FORMAT (1H  ) 

6 FORMAT (15. 5F10.3) 

7 FORMAT (1H1) 

9  FORMAT ( 515, 2F5. 2, 15) 

10  FORMAT (55H  I EL  WIDTH-FT  HEIGHT-FT  LENGTH-FT  OTHER  SELECT  ) 
1 FORMAT (A80) 

READ (5,1)  TITLE 

READ (5. 6)  (NEXT ( IC) .WDTH(IC) , HGHT(IC) .XLTH(IC) .OTHR(IC) , 

* SELECT ( IC ) ,IC=1. NEL) 

IF ( IPRINT . EQ . 0 ) GO  TO  500 
WRITE (6, 5) 

WRITE (6,1)  TITLE 
WRITE (6. 5) 

WRITE (6. 10) 

WRITE (6. 6) ( NEXT ( IC ) .WDTH(IC) , HGHT(IC) .XLTH(IC) .OTHR(IC) . 

* SELECT ( IC ) ,IC=1. NEL) 

500  CONTINUE 

COF ( 1 ) =WDTH ( 1 ) 

PEN=HGHT ( 1 ) 

ZDOOR=XLTH ( 1 ) 

53  FORMAT ( 55H  l=DOOR,  PEN  WALL  (FT),  DOOR  THICKNESS (FT  OF  CONCRETE)  ) 
IF ( IPRINT . LT . 4 ) GO  TO  501 
WRITE (6. 5) 

WRITE (6, 53) 

WRITE (6, 6)  NEXT(l) .PEN. ZDOOR 
WRITE (6, 5) 

501  CONTINUE 
ECC ( 1 ) =PEN 
DELR(l) = ZDOOR 
AREA ( 1 ) = 0 . 

DO  25  IC=2 , NEL 
IEL=NEXT ( IC) 

IF ( IEL-7)  11,30,33 

11  AREA ( IC) =WDTH ( IC) *HGHT ( IC ) 

RAT= (2. *PEN+WDTH(IC) ) * ( 2 . *PEN+HGHT ( IC) ) /AREA(IC) 

SCL ( IC ) =SQRT ( RAT) 

ECC(IC) =WDTH(IC) /HGHT(IC) 

DZRO ( IC) =OTHR ( IC ) 

DELR( IC) =1 . 

IF ( IEL . EQ . 3 ) GO  TO  13 
DELR(IC) =OTHR(IC) 

IF ( OTHR ( IC ) . EQ . 0 . ) DELR(IC)=1. 

IF ( IEL . EQ . 6 ) GO  TO  22 
13  COF ( IC )=XLTH(IC)/SQRT( AREA (IC) ) 

GO  TO  25 

22  COF ( IC )= WDTH (IC) 

ECC ( IC ) =HGHT ( IC ) 

DELR ( IC) =XLTH ( IC ) 

DZRO ( IC ) =OTHR ( IC ) 

SCL(IC) =OTHR ( IC+1 ) 

GO  TO  25 
30  ICP=IC-1 
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277 
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286 
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289 

290 

291 

292 

293 

294 

295 

296 

297 
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299 

300 

301 

302 

303 

304 
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314 

315 
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322 
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325 


FILE  1 


IF (NEXT ( IC-l ) . EQ . 8 ) ICP=IC-3 
COF ( IC) =XLTH ( IC ) /SQRT (AREA ( ICP ) ) 

SCL ( IC ) = AREA ( 2 ) /AREA ( ICP ) 

IF (NEXT ( IC-l ) . EQ . 8 ) SCL ( IC ) =SCL ( IC-l ) 

DELR ( IC) =DELR ( ICP ) 

ECC ( IC) =ECC ( ICP) 

DZRO ( IC) =XLTH ( IC) - . 01 
GO  TO  25 

33  COF ( IC ) = WDTH ( IC ) 

ECC(IC)=HGHT(IC) 

DELR ( IC ) =XLTH ( IC ) 

DZRO ( IC ) =OTHR ( IC ) /AREA ( IC- 2 ) 

SCL ( IC ) = AREA ( 2 ) /OTHR ( IC ) 

DZRO ( 1 ) =OTHR ( 1 ) 

25  CONTINUE 
1000  RETURN 
END 

SUBROUTINE  SRSDET ( ISORS . ING . MODG , MODN ) 

COMMON  DIAMS ( 18 ) . DUCTO ( 3 . 17 . 6 ) .00011(3,17,6) , DUCT2 (3,17,6) .IEL.LP 
X BEND (10,6) , SGRUPS ( 23 ) , DWTS ( 21 ) , CRRNTO (9,10) , RCRRNT (9,10) 
COMMON  RVCT ( 21 ) . RMTX ( 21 . 21 ) , SO ( 9 ) , SP ( 9 ) , DZRO ( 9 ) , NEXT ( 9 ) , COF ( 9 ) , 

X ECC ( 9 ) , DELR ( 9 ) ,SNYNEU(23) ,SINVAL(10) , C01 , C02 , VEL (9,21) , IPRINT 
COMMON  RCOSO , RCOSP , RCOS , DGRUPS ( 23 ) , VELO (9,21). VELP (9,21) 

X ,VEL1(5) .SCALE. SWTS(23. 21) ,WTSD(23,21) 

COMMON  IDOSE . HNDRSN ( 23 ) , SEPI ( 23 ) . STH ( 23 ) , SINR ( 23 ) , THMLN ( 23 ) 

X , ENGAM ( 21 ) , DIMEN , SCOBLT ( 23 ) . AWTS (6.5), SGAM ( 23 ) , DGAM ( 23 ) , AMOD (5,5) 
DIMENSION  SWTD ( 21 ) , TEMP ( 5 ) 

1 FORMAT (1H  ) 

2 FORMAT ( 5 5H  DWTS //SWTD/ /ENGAM 

3 FORMAT ( 7 ( 1PE10 . 3 ) ) 

CNG=ING 

CG=CNG**2 

CN=1 . + (CNG-CG) /2 . 

CN=CN*DIMEN 
CG=CG* DIMEN 

JSTT=0 

IOX=0 

CALL  WTLISTS ( JSTT . JXWT . IOX . NMAX . DWTS . SWTD , ENGAM , IDOSE . ISORS , LP ) 

IF ( IPRINT, GT. 2) GO  TO  100 
IF ( JSTT . EQ . 0 ) GO  TO  15 

50  IF ( IDOSE.LT. 5)  GO  TO  100 
IAG=IS0RS-3 

IF(IAG.LT.l)  GO  TO  100 
XL=COF ( 2) 

CUT=0 . 5/SQRT (XL**2+0 , 25 ) 

CUTP=SQRT ( 1 . -CUT**2 ) 

DO  97  J=1 , 5 
AMOD (3, J)=0. 

97  AMOD ( 4 , J) =0 . 

DO  98  J=1 , 5 
FJ= J 

IF ( CUT . GT . AMOD ( 1 , J ) ) GO  TO  99 

98  CONTINUE 

99  NCUT=FJ+ . 1 
DO  980  J=1 . 5 
FJ=  J 

IF ( CUTP  GT . AMOD ( 1 , J ) ) GO  TO  990 
980  CONTINUE 
990  NCUTP=FJ+ . 1 

FR= ( AMOD ( 2 , NCUT ) -CUT ) / ( AMOD ( 2 . NCUT ) -AMOD ( 1 , NCUT ) ) 

FRP= ( AMOD ( 2 . NCUTP ) -CUTP ) / ( AMOD ( 2 , NCUTP ) -AMOD ( 1 . NCUTP ) ) 

DO  96  J=1 ,NCUT 
96  AMOD ( 3 , J) =AWTS ( IAG. J) 

DO  960  J=l, NCUTP 
960  AMOD ( 4 , J ) = AWTS ( IAG+  2 . J ) 

AMOD ( 3 , NCUT ) =FR*AMOD ( 3 . NCUT ) 

AMOD ( 4 . NCUTP ) =FRP* AMOD ( 4 . NCUTP ) 

SD  = 0 . 

ED=0 . 

DO  102  J=1 , 5 
SD=SD+AMOD ( 3 . J) 

102  ED=ED+AMOD (4 . J) 

DO  103  J=1 , 5 

AMOD ( 3 . J) = AMOD (3,J)/(SD+l.E-7) 

103  AMOD ( 4 , J) =AMOD ( 4 , J) / ( ED+1 . E-7 ) 

TOT=0 . 

TVP=VELP (2,2) + VELP (2.3) + VELP (2,4) +VELP (2.5) +VELP (2,6) +VELP (2.7) 
TVO=VELO (2,1) 

IF ( IEL . LT . 4 ) TVO=0 . 

DO  200  J=1 , 5 
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FILE  1 


Continued 


TEMP ( J ) =TVO*AMOD ( 4 . J) +TVP*AMOD ( 3 . J) 

200  TOT  = TOT + TEMP ( J ) 

DO  201  J=1 , 5 

201  VEL1( J) =TEMP(6-J) /TOT 

IF (IPRINT.LT. 4)  GO  TO  100 
WRITE (6,1) 

WRITE (6. 2) 

WRITE (6,3)  ( AMOD (3,J) ,J=1,5) 

WRITE (6,3)  ( AMOD (4,J) ,J=1,5) 

WRITE (6. 3)  TVO , TVP , TOT . TOT , (VEL1(J) , J=l,5) 

WRITE (6, 2) 

498  WRITE (6,1) 

100  DO  12  MEG=1 , 23 

GO  TO  (105,110,115,120,125) .ISORS 
105  SGRUPS (MEG) =SINR (MEG) 

GO  TO  4 

110  SGRUPS (MEG) =STH( MEG) 

GO  TO  4 

115  SGRUPS (MEG) =SEPI (MEG) 

GO  TO  4 

120  SGRUPS (MEG) =SGAM( MEG) 

GO  TO  4 

125  SGRUPS ( MEG )=SCOBLT( MEG) 

4 GO  TO  (205.210,215,220,225) , IDOSE 
205  DGRUPS ( MEG ) = SNYNEU ( MEG ) 

GO  TO  12 

210  DGRUPS (MEG) =1. 

GO  TO  12 

215  DGRUPS (MEG) =HNDRSN( MEG) 

GO  TO  12 

220  DGRUPS ( MEG ) =THMLN( MEG) 

GO  TO  12 

225  DGRUPS (MEG) =DGAM( MEG) 

12  CONTINUE 

IF ( IDOSE . LT . 5 ) GO  TO  13 

NMAX=10 

MMAX=17 

DO  14  N=1 , 21 

DWTS (N) =0 . 

SWTD ( N) =0 . 

14  ENGAM (N) =0 . 

CALL  SGD ( SGRUPS , SWTS , WTSD , DGRUPS , VEL1 , MODG ) 

GO  TO  11 

13  CALL  SND( SGRUPS, SWTS, WTSD. DGRUPS, ENGAM, MODN) 

NMAX=21 

MMAX=23 

11  IF ( JSTT . EQ . 0 ) GO  TO  15 
DO  19  N=1,NMAX 
DWTS (N) =0 . 

SWTD ( N ) = 0 . 

DO  19  M=1 , MMAX 

SWTD(N) =SWTD (N) +SWTS(M,N) * DGRUPS (M) 

DWTS ( N ) =DWTS ( N ) + WTSD ( M . N ) * SGRUPS ( M ) 

19  CONTINUE 

15  CALL  WTLISTS ( JSTT, JXWT, IOX. NMAX. DWTS, SWTD, ENGAM, IDOSE, ISORS, LP) 
DO  20  N=1 , NMAX 

IF ( IDOSE. EQ. 5)  DWTS (N) =SWTD (N) 

DWTS  (N)  = CN* DWTS  (N)  + CG* ENGAM ( N ) 

20  CONTINUE 
C 

25  IF (IPRINT.LT. 2)  GO  TO  1000 
WRITE (6.1) 

WRITE (6. 2) 

WRITE (6,3)  (DWTS(N) ,N=1.NMAX) 

WRITE (6.1) 

WRITE (6,3)  (SWTD(N) ,N=1,NMAX) 

WRITE (6.1) 

WRITE (6,3)  ( ENGAM (N) ,N=1,NMAX) 

WRITE (6. 1) 

1000  RETURN 
END 

SUBROUTINE  SND ( SGRUPS , SWTS . WTSD . DGRUPS . ENGAM . MODN ) 

SLARGE  ALBDO 

DIMENSION  SGRUPS (23) .SWTS (23. 21) .WTSD (23. 21) .DGRUPS (23) . 

* ALBDO (23. 21. 23) .ENGAM (21) 

CHARACTER* 80  TITLE 
1 FORMAT (A80) 

7010  FORMAT ( 15 ) 

7020  FORMAT ( 8 ( 1PE10 . 3 ) ) 

IF  ( MODN . EQ  . 1 ) GO  TO  12 

IUT=11 

REWIND  IUT 

READ ( IUT . 1 ) TITLE 
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IF (IPRINT.LT. 4)  GO  TO  12 
WRITE (6.1)  TITLE 
12  FG=. 095/. 75*4. 5E-10 
DO  25  N=1 , 21 
ENGAM (N) =0 . 

DO  25  M=1 » 23 
SWTS(M.N) =0. 

25  WTSD (M , N) =0 . 

DO  350  MEG=1 , 23 
IF(MODN.EQ.l)  GO  TO  35 
DO  5 M=1 , 20 
DO  5 J=1 , 23 
5 ALBDO ( J , M , MEG ) = 0 . 

3 READ (IUT. 7010)  ISORS 
DO  30  J=ISORS . 23 

30  READ (IUT, 7020)  (ALBDO (J.M. MEG) , M=l, 20) 

35  SWTS ( ISORS, 1 ) =1 . *SGRUPS ( ISORS) 

WTSD ( ISORS . 1 ) =1 . *DGRUPS ( ISORS ) 

DO  100  N=2 , 21 
DO  100  M=ISORS , 23 

SWTS (M . N) =SWTS (M . N) + ALBDO (M.N-l, MEG) *SGRUPS ( ISORS ) 

ENGAM ( N ) = ENGAM ( N ) + ALBDO (M.N-l. MEG ) *FG*  SGRUPS ( ISORS ) 

100  WTSD ( ISORS. N)=WTSD( ISORS. N)  + ALBDO (M.N-l .MEG) *DGRUPS (M) 
350  CONTINUE 

IF(MODN.EQ.O)  CLOSE(ll) 

MODN=l 

RETURN 

END 

SUBROUTINE  SGD ( SGRUPS , SWTS . WTSD , DGRUPS , VEL1 , MODG ) 

$ LARGE  ALB 

DIMENSION  SGRUPS ( 23 ) , SWTS (23,21), WTSD (23,21). DGRUPS ( 23 ) , 
* ALB (17,9,85) ,VEL1(5) 

CHARACTER* 80  TITLE 

1 FORMAT (A80) 

2 FORMAT (215) 

3 FORMAT ( IX. 7 ( 1PE9 . 3 ) , 2 ( 1PE8 .2)) 

IF (MODG. EQ . 1 ) GO  TO  12 
IUT=10 

REWIND  IUT 
READ ( IUT. 1)  TITLE 
IF (IPRINT.LT. 4) GO  TO  12 
WRITE (6,1)  TITLE 
12  CONTINUE 

DO  25  M=1 , 23 
DO  25  N=1 , 21 
SWTS(M.N) =0. 

25  WTSD (M , N) =0 . 

DO  350  K=1 , 85 
IF (MODG . EQ . 1 ) GO  TO  35 
KE=(K-l)/5 
LOB=K-5*KE 
DO  5 M=1 . 9 
DO  5 J=1 . 17 
5 ALB ( J , M , K) =0 . 

READ ( IUT , 2 ) IK. IE 
DO  30  J=IE. 17 

30  READ (IUT, 3)  ( ALB ( J. M , K) , M=1 . 9 ) 

35  SWTS (IE. 1)  =SWTS ( IE , 1 ) +VEL1 ( LOB ) *SGRUPS ( IE ) 

WTSD ( IE . 1 ) = WTSD ( IE . 1 ) +VEL1 ( LOB ) *DGRUPS ( IE ) 

DO  100  N=2 , 10 
DO  100  M=IE , 17 
IF ( N . LE . 2 ) GO  TO  40 
IF ( M . GT . 16 ) GO  TO  40 
IF(N.GT.M-IO)  GO  TO  100 
40  CONTINUE 

SWTS(M.N)  =SWTS  (M , N)  + ALB  (M . N-l , K)  * SGRUPS  ( IE)  *VEL1  (LOB ) 
100  WTSD ( IE . N ) =WTSD ( IE , N) +ALB (M . N-l , K ) *DGRUPS (M ) *VEL1 ( LOB ) 
350  CONTINUE 

IF (MODG . EQ . 0 ) CLOSE (10) 

MODG=l 

RETURN 

END 

SUBROUTINE  ORDST ( RHO , A . AL , ORLIST ) 

DIMENSION  TEMP (21,21). BLST ( 21 ) , ORLIST ( 21 ) 

AB=A*RHO 

B= ( A-AB ) / ( 1 . -AB* (1. 36788-. 36788*A) ) 

ALB=  ( 1 . -B ) * ( 1 . -AB ) *AL/RHO/ ( 1 . -A) 

BLST ( 1 ) =1 . 

TEMP (1.1)  =1. 

DO  10  1=2.21 

BLST ( I ) =B*BLST ( 1-1 ) 

10  TEMP (1,1) =TEMP ( 1 , 1-1 ) *AB 
TEMP ( 1 , 1 ) =0 . 
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DO  20  J=2 , 21 
DO  20  I=J, 21 

20  TEMP ( J, I ) =TEMP ( J-l , 1-1 ) 

DO  30  1=1.21 
ORLIST ( I ) =0 . 

DO  30  J-l.I 

30  ORLIST ( I ) =ORLIST ( I ) + TEMP ( J, I) *BLST ( J ) 

DO  40  1=1.21 

40  ORLIST ( I ) =ORLIST ( I ) *ALB 
RETURN 
END 

SUBROUTINE  DCTISO ( DLNGTH , IC ) 

COMMON  DIAMS ( 18 ) , DUCTO ( 3 , 17 . 6 ) . DUCT1 ( 3 . 17 , 6 ) , DUCT2 (3.17,6) .IEL.LP, 
X BEND (10,6), SGRUPS ( 23 ) , DWTS ( 21 ) , CRRNTO (9,10), RCRRNT (9,10) 
COMMON  RVCT ( 21 ) , RMTX ( 21 , 21 ) . SO ( 9 ) , SP ( 9 ) , DZRO ( 9 ) , NEXT ( 9 ) , COF ( 9 ) , 

X ECC ( 9 ) , DELR ( 9 ) ,SNYNEU(23) .SINVAL(IO) , COl , C02 , VEL (9,21) , IPRINT 
COMMON  RCOSO , RCOSP , RCOS , DGRUPS ( 23 ) , VELO (9,21), VELP (9.21) 

X , VEL1 ( 5 ) . SCALE , SWTS (23,21). WTSD (23,21) 

COMMON  I DOSE , HNDRSN ( 23 ) , SEPI ( 23 ) , STH ( 23 ) , SINR ( 23 ) , THMLN ( 23 ) 

X , ENGAM ( 21 ) , DIMEN , SCOBLT ( 23 ) , AWTS (6,5), SGAM ( 23 ) . DGAM ( 23 ) . AMOD (5,5) 
DIMENSION  XI ( 5 ) , ORLIST ( 21 ) , ESC ( 21 ) 

DO  4 1=1.5 
4 XI (I) =0. 

1 FORMAT (1H  ) 

2 FORMAT ( 13 , 11 ( 1PE10 . 2 ) ) 

3 FORMAT ( 13 , 7E13 . 3 ) 

IF ( DELR ( IC ) . LT . 0 . ) GO  TO  131 
DO  130  J=1 , 21 

130  ESC( J) =DELR(IC) 

GO  TO  133 

131  ESC(l) =1 . 

DO  132  J=2 , 21 

132  ESC ( J) =ESC ( J-l ) * ( -DELR ( IC ) ) 

133  CONTINUE 

DLNGTH =DLNGTH /SCALE 
EC=ECC(IC) 

ECCO= - . 06  * LOG ( EC ) * * 2 * DLNGTH / SQRT (3.182  + DLNGTH  * * 2 ) 

ECCP=- . 109*L0G(EC) **2 
ECC0=EXP ( ECCO ) 

ECCP=EXP ( ECCP ) 

XI ( 1 ) =2 . 

XI  ( 4 ) =2 . 

IF  ( DLNGTH . GE . 8 . ) GO  TO  10 
LX= ( 2 . *DLNGTH+1 . ) 

PWR= (DLNGTH-DIAMS(LX) )/( DIAMS ( LX+1 ) -DIAMS ( LX)  ) 

DO  15  KX=1 . 5 

Y2= DUCTO (LP , LX+1 , KX) 

Yl= DUCTO ( LP , LX, KX) 

RCRRNT (IC.KX) =0. 

15  CRRNTO ( IC. KX) =Y1* (Y2/Y1) **PWR 
GO  TO  13 
10  DO  12  KX=1 , 5 

AO=DUCTO (LP, 17 , KX) 

A2= DUCTO ( LP , 13 , KX) 

CALL  XTRAP ( AO . A2 , XI ( KX ) , DLNGTH , AXTR ) 

12  CRRNTO (IC,KX)= AXTR 

13  CONTINUE 

AL= ( - . 35+CRRNTO ( IC, 2 ) )/( . 5+DLNGTH**3 ) 

A=- . 30+CRRNTO ( IC , 3 ) 

RHO= . 95 

RMTX (1,1) = CRRNTO ( IC , 1 ) 

CALL  0RDST(RH0. A, AL, ORLIST) 

DO  19  1=2,21 
RMTX (1,1) =ORLIST ( I ) 

DO  19  J=2 , I 

19  RMTX ( J , I ) =RMTX ( J-l , 1-1 ) 

20  SO ( IC ) = 1 . -CRRNTO ( IC , 4 ) 

SP (IC) =-l . +2 . *CRRNTO ( IC , 5 ) 

VELO ( IC , 1 ) =RMTX (1,1) *RVCT ( 1 ) *ESC ( 1 ) *ECCO 
DO  31  J=2 , 21 

VELO ( IC , J ) =RMTX ( J, J) *RVCT ( J ) *ESC ( 1 ) *ECCO 
VELP ( IC , J) =0 . 

IP=J-1 

DO  30  1=1, IP 

30  VELP ( IC , J) =VELP ( IC . J) +RMTX ( I , J) *RVCT ( I ) *ESC ( J-I+l ) 

VELP ( IC , J ) =VELP ( IC , J) *ECCP 

31  CONTINUE 

DO  32  J=1 . 21 

VEL ( IC , J ) = VELO ( IC , J ) +VELP ( IC, J) 

32  RVCT ( J ) = VEL ( IC . J ) 

DLNGTH  = DLNGTH  * SCALE 

IF (IPRINT.LT. 4)  GO  TO  1000 
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40  WRITE (6.1) 

WRITE (6,3)  IC  , DLNGTH , ( CRRNTO (IC.KX) ,KX=1,5) 
WRITE (6,3)  IC.  CRRNTO ( IC , 1 ) , A , AL . SO ( IC ) , SP ( IC ) 
WRITE  (6, 2)  I EL,  (VELOUC.  J)  . J=l.ll) 

IND=1 

WRITE (6, 2)  IND. (VELP(IC. J) , J=l.ll) 

IND  = 2 

WRITE (6,2)  IND, (RVCT( J) , J=l,ll) 

1000  RETURN 
END 
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SUBROUTINE  DCTSLT ( DLNGTH . IC ) 

COMMON  DIAMS ( 18 ) . DUCTO ( 3 . 17 , 6 ) , DUCT1 ( 3 . 17 . 6 ) , DUCT2 (3,17,6) .IEL.LP. 
X BEND (10.6). SGRUPS ( 23 ) . DWTS ( 21 ) , CRRNTO (9.10). RCRRNT (9.10) 

COMMON  RVCT ( 21 ) . RMTX ( 21 . 21 ) . SO ( 9 ) . SP ( 9 ) . DZRO ( 9 ) . NEXT ( 9 ) . COF ( 9 ) , 

X ECC ( 9 ) , DELR ( 9 ) . SNYNEU ( 23 ) , SINVAL ( 10 ) , COl . C02 . VEL (9.21). 1PRINT 
COMMON  RCOSO . RCOSP . RCOS . DGRUPS ( 23 ) , VELO (9.21). VELP (9,21) 

X , VEL1 ( 5 ) . SCALE , SWTS (23.21). WTSD (23.21) 

COMMON  IDOSE , HNDRSN ( 23 ) , SEPI ( 23 ) , STH ( 23 ) . SINR ( 23 ) , THMLN ( 23 ) 

X , ENGAM ( 21 ) , DIMEN , SCOBLT ( 23 ) , AWTS(6.5) ,SGAM(23) ,DGAM(23) ,AMOD(5,5) 
DIMENSION  RMTX1 (21,21) , RMTX2 ( 21 . 21 ) ,V01(21) ,V02(21) ,VP1(21) . 

X VP2 ( 21 ) . XI ( 5 ) . ORLIST ( 21 ) , ESC ( 21 ) 

DO  4 1=1,5 
4 XI (I) =0. 

1 FORMAT (1H  ) 

2 FORMAT ( 13 , 11 ( 1PE10 .2)) 

3 FORMAT ( 13 . 9E12 . 4 ) 

IF ( DELR ( IC ) . LT . 0 . ) GO  TO  131 
DO  130  J=1 , 21 

130  ESC( J) =DELR(IC) 

GO  TO  133 

131  ESC ( 1 ) =1 . 

DO  132  J=2 , 21 

132  ESC ( J) =ESC ( J-l ) * ( -DELR ( IC) ) 

133  CONTINUE 

DLNGTH  = DLNGTH /SCALE 
C01= . 75/ (1 . +SO(IC-l) ) -0 . 5 
C02=l . -COl 
U= . 5* ( 1 . +SO ( IC-1 ) ) 

ZM=U/SQRT ( 1 . -U**2 ) 

VALO=l . - DLNGTH /ZM 
IF ( VALO . LT . 0 . ) VALO=0 . 

7 EC=ECC ( IC) 

ECCO  = -.06*  LOG (EC) **2* DLNGTH / SQRT (3.182  + DLNGTH* * 2 ) 

ECCP=- . 109*LOG(EC) **2 
ECCO=EXP(ECCO) 

ECCP=EXP ( ECCP) 

IF ( DLNGTH. GE. 8. ) GO  TO  10 
LX= ( 2 . *DLNGTH+ 1 . ) 

PWR= (DLNGTH-DIAMS (LX) )/( DIAMS ( LX+1 ) -DIAMS ( LX) ) 

DO  15  KX=1 , 5 
Y2=DUCT1 (LP.LX+l.KX) 

Y1=DUCT1 ( LP . LX , KX) 

RCRRNT ( KX . 1 ) =Y1 *( Y2 /Y1 )** PWR 
Y2=DUCT2( LP.LX+l.KX) 

Y1=DUCT2 ( LP , LX , KX) 

15  RCRRNT  ( KX,  2 ) = Y1*  ( Y2/Y1 ) “PWR 
GO  TO  13 
10  DO  12  KX=1 , 5 

A0=DUCT1 ( LP , 17 , KX ) 

A2=DUCT1 ( LP , 13 , KX) 

CALL  XTRAP (AO.A2,XI(KX), DLNGTH . AXTR ) 

RCRRNT ( KX. 1 ) =AXTR 
AO=DUCT2(LP. 17, KX) 

A2=DUCT2 ( LP . 13 , KX) 

CALL  XTRAP ( AO , A2 . XI ( KX ) . DLNGTH , AXTR ) 

12  RCRRNT (KX. 2) =AXTR 

13  CONTINUE 

AL1= ( - . 35+RCRRNT ( 2 , 1 ) ) / ( . 5+DLNGTH**3 ) 

AL2= ( - . 35+RCRRNT (2,2) ) /( . 5+DLNGTH**3 ) 

Al=- . 3+RCRRNT (3.1) 

A2=- . 3+ RCRRNT (3,2) 

RHO= . 85 

RMTX1 (1.1) =VALO 
RMTX2 (1,1) =VALO 
CALL  0RDST(RH0.A1.AL1. ORLIST) 

DO  17  1=2.21 
17  RMTX1 (1,1) =ORLIST ( I ) 

CALL  ORDST(RHO,A2.AL2. ORLIST) 

DO  19  1=2.21 
RMTX2 (1,1) =ORLIST ( I ) 

DO  19  J=2 , I 

RMTX1( J.I) =RMTX1( J-l.I-l) 

19  RMTX2 (J.I) =RMTX2 ( J-l , 1-1 ) 

20  RC01= 2. /RCRRNT (4.1) 

RC02=2 . / ( 1 . +RCRRNT (4,2) ) 

RCP1=1. /RCRRNT (5,1) 

RCP2=1 . /RCRRNT (5.2) 

V01T=1 . 010 IE- 9 
V02T=2 . 002E-9 
VP1T=1 . 0101E-9 
VP2T=2 . 0202E-9 
DO  31  J=1 , 21 
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137: 
138: 
139: 
140: 
141: 
142: 
143: 
144  : 
145: 
146: 
147: 
148: 
149: 
150: 
151: 
152: 
153: 
154; 
155: 
156: 
157: 
158: 
159: 
160: 
161: 
162: 
163: 
164: 
165: 


V01 ( J ) =RMTX1 ( J . J) *VELO ( IC-1 . J ) *ESC ( 1 ) 
V02 ( J) =RMTX2 ( J,  J) *VELO ( IC-1 , J) *ESC ( 1) 
VP1 ( J) =0. 

VP2 ( J) =0 . 

JP=J-1 


IF(J-l)  28,28,29 

29  DO  30  1=1, JP 

VP1(J)=VP1( J)+RMTX1(I, J)*VELO(IC-l.I)*ESC( J-Id) 

30  VP 2 ( J) =VP2 ( J) +RMTX2 ( I , J) *VELO ( IC-1 , I ) *ESC ( J-I+l ) 

28  V01T=V01T+V01 ( J) 

V02T=V02T+V02 ( J) 

VP1T=VP1T+VP1 ( J) 

VP2T=VP2T+VP2 ( J) 

31  CONTINUE 

DO  32  J=1 , 21 

VELO ( IC, J) =C01*V01 ( J) +C02*V02 ( J) 

32  VELP ( IC , J) =C01*VP1 ( J) +C02*VP2 ( J) 

DO  33  J=1 , 21 

VELO(IC, J) =VELO(IC, J) *ECC0 
VELP ( IC . J ) = VELP ( IC , J ) *ECCP 
VEL ( IC , J) =VELO ( IC , J) +VELP ( IC , J ) 

33  RVCT ( J ) = VEL ( IC , J ) 

SO(IC) =-l. +2.* (C01*V01T+C02*V02T)/(C01*V01T*RC01+C02*V02T*RC02) 

SP ( IC) =-l . +2 . * (C01*VP1T+C02*VP2T) / (C01*VP1T*RCP1+C02*VP2T*RCP2 ) 
DLNGTH  = DLNGTH  * SCALE 
IF ( IPRINT . NE . 4 ) GO  TO  1000 
40  WRITE (6,1) 

IC  .DLNGTH, (RCRRNT(KX.l)  ,KX=1,5) 

IC  .DLNGTH. (RCRRNT(KX, 2)  ,KX=1,5) 

IC , AL1 , A1 , AL2 , A2 , COl , SO ( IC ) .SP(IC) ,C02 
IC . V01T . V02T , VP1T , VP2T . RCOl . RC02 . RMTX1 (1.1), RMTX2 (1.1) 
IEL, (VELO(IC. J) , J=l,ll) 

IEL, (VELP ( IC. J) . J=l.ll) 


IND. (VOl(J) , J=l,21.2) 
IND. ( V02 ( J) , J=l,21.2) 


WRITE (6.3) 

WRITE (6. 3) 

WRITE (6. 3) 

WRITE (6, 3) 

WRITE (6. 2) 

WRITE (6, 2) 

IND=13 
WRITE (6. 2) 

WRITE (6, 2) 

1000  RETURN 
END 

SUBROUTINE  XTRAP ( AO . A2 . XI , FA . AINT ) 

SAVE=XI 

FO  = 8 . 

F2  = 6 . 

XI=L0G(A2/A0) /LOG ( FO/F2 ) - . 05 

DL=0 . 1 

FOL=FO**DL 

F2L=F2**DL 

FAL=FA**DL 

DO=AO*FO**XI*FOL 

D2=A2*F2**XI*F2L 


AL1= ( DO-D2 ) / ( FOL-F2L ) 

IF(ALl)  10.10.5 

5 AL2= ( F0L*D2-F2L*D0 ) / (FOL-F2L) 

AINT= ( AL1+AL2/FAL) /FA**XI 
GO  TO  1000 

10  XI=LOG ( A2/AO) /LOG ( FO/F2 ) 

ALX= ( AO*FO** (XI+1 . ) -A2*F2** (XI+1 ) )/(FO-F2) 

AINT=AL1/FA**XI 
XI = SAVE 
1000  RETURN 
END 

SUBROUTINE  TURN ( SEN. IC) 

COMMON  DIAMS (18 ) . DUCTO (3.17.6), DUCT1 (3.17,6), DUCT2 (3.17,6), IEL . LP . 
X BEND (10.6), SGRUPS ( 23 ) , DWTS ( 21 ) , CRRNTO (9.10). RCRRNT (9.10) 

COMMON  RVCT ( 21 ) , RMTX ( 21 . 21 ) , SO ( 9 ) , SP ( 9 ) . DZRO ( 9 ) . NEXT ( 9 ) , COF ( 9 ) . 

X ECC ( 9 ) . DELR ( 9 ) . SNYNEU ( 23 ) . SINVAL ( 10 ) . COl , C02 . VEL (9.21), IPRINT 
COMMON  RCOSO , RCOSP . RCOS , DGRUPS ( 23 ) . VELO (9,21), VELP (9,21) 

X , VEL1 ( 5 ) , SCALE . SWTS (23.21), WTSD (23.21) 

COMMON  I DOSE , HNDRSN ( 23 ) . SEPI ( 23 ) , STH ( 23 ) . SINR ( 23 ) , THMLN ( 23 ) 

X , ENGAM ( 21 ) . DIMEN , SCOBLT ( 23 ) , AWTS (6,5), SGAM ( 23 ) . DGAM ( 23 ) , AMOD (5.5) 
DIMENSION  RX1 ( 13 ) , RINF (13 ) , ORLIST ( 21 ) , ESC ( 21 ) 

DATA  RX1/.88, .86, .84, .82,  79, . 74 , . 69 , . 62 , . 54 . .44, .26, .01, .001/ 
DATA  RINF/ . 86 . . 83 . . 80 , . 76 . . 71 , . 66 , . 58 , . 48 , . 37 , . 24 , . 01 , 

* .001.. 0001/ 

IF ( DELR (IC ) . LT . 0 . ) GO  TO  131 
DO  130  J=1 , 21 

130  ESC ( J) =DELR (IC ) 

GO  TO  133 

131  ESC (1 ) =1 . 

DO  132  J=2 , 21 

132  ESC ( J) =ESC ( J— 1 ) * ( -DELR (IC) ) 

133  CONTINUE 


C0=0. 


CP=0. 
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172 
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185 
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197 
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199 
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201 
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203 
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207 
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210 
211 
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213 

214 

215 

216 

217 

218 
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222 

223 
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236 
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240 
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244 

245 
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AZFAC=1 . -DZRO ( IC ) 

CP1=0. 

CP2=0 . 

RC0S0=2 . / ( 1 . +S0 ( IC-1 ) ) 

RC0SP=2 . / ( 1 . +SP ( IC-1) ) 

DO  305  1=1.21 
CO=CO+VELO ( IC-1 , I ) 

CP=CP+VELP ( IC-1 . I ) 

VELO ( IC , I ) =0 . 

VELP ( IC , I ) =0 . 

305  RVCT ( I ) = VELO ( IC-1 . I ) 

RCOSO= ( CO*RCOSO+CP*RCOSP ) / ( CO+CP ) 

RCOSP=RCOSO 

RCOS=RCOSO 

IX=1 

134  SEN=2 . / (RCOS+1 . E-6 ) -1 . 

C 

RCSS=COF (IC) 

KS= ( SEN* 10 . +1 . ) 

IF(KS.LT.l)  KS=1 
RKS=KS 

PWR= ( SEN*10 . +1 . 00001-RKS ) 

C1=RX1 ( KS ) * (RX1 ( KS+1 ) /RX1 ( KS ) )**PWR 
CINF=RINF ( KS ) * ( RINF (KS  + 1) /RINF ( KS ) ) * *PWR 
BD= ( 1 . -Cl ) / (Cl-CINF ) 

BN=CINF*BD 
RX=RCSS**1 . 5 

RDUX= ( 1 . +BN*RX) / ( 1 . +BD*RX) 

C 

1 FORMAT (1H  ) 

2 FORMAT ( 13 , 11 ( 1PE10 . 2 ) ) 

3 FORMAT ( 13 , 7E13 . 3 ) 

LS= ( 5 . *SEN+3 . ) 

DO  15  KX=1 , 5 
Y2=BEND ( LS+1 , KX ) 

Y1=BEND ( LS . KX) 

PWR= ( SEN-SINVAL ( LS ))/( SINVAL ( LS +1 ) -SINVAL ( LS )+. 00001 ) 
15  CRRNTO ( IC, KX) =Y1* ( Y2/Y1) **PWR 
AL=CRRNTO ( IC , 2 ) 

A=CRRNTO ( IC , 3 ) 

SO ( IC) =CRRNTO ( IC , 4 ) *SQRT ( 2 . -SEN) -1 . 

SP ( IC) =2 . * CRRNTO ( IC . 5 ) -1 . 

IF(IX.GT.O)  S01=S0(IC) 

IF(IX.GT.O)  SP1=SP(IC) 

IF(IX.LT.O)  S02=S0 ( IC) 

IF(IX.LT.O)  SP2=SP(IC) 

RHO= . 75 

RMTX (1.1) = CRRNTO ( IC . 1 ) 

CALL  ORDST ( RHO , A , AL , ORLIST ) 

DO  20  1=2.21 

RMTX ( 1 . I ) =ORLIST ( I ) 

DO  19  J=2 . I 

19  RMTX ( J, I ) =RMTX ( J-l . 1-1 ) 

20  CONTINUE 

VELO ( IC . 1 ) =VELO ( IC , 1 ) +RMTX (1.1) *RVCT ( 1 ) *ESC ( 1 ) 
IF(IX.GT.O)  C01=VEL0 ( IC , 1 ) 

IF(IX.LT.O)  C02=VEL0 ( IC, 1 ) 

DO  31  J=2 , 21 

VELO(IC. J) =VELO(IC. J) +RMTX( J, J) *RVCT ( J ) *ESC ( 1 ) 

IP=J-1 

DO  30  1=1, IP 

DVP=RMTX ( I . J) *RVCT ( I ) *ESC ( J-I+l ) *RDUX 

30  VELP ( IC , J) =VELP ( IC , J) +DVP 
IF(IX.GT.O)  C01=C01+VEL0 ( IC , J) 

IF(IX.GT.O)  CPI =CP1+ VELP ( IC , J) 

IF(IX.LT.O)  C02=C02+VEL0 ( IC , J) 

IF(IX.LT.O)  CP2=CP2+VELP(IC, J) 

31  CONTINUE 
C 

IF(IX.LT.O)  GO  TO  325 

REN0RM=1 . + ( 1 . -AZFAC) *C01/ (C0+CP*1 . E-5 ) 

DO  315  J=1 , 21 

VELO(IC, J) =VELO ( IC , J) * AZFAC 
VELP ( IC . J ) = VELP ( IC , J ) *RENORM 
315  RVCT ( J) =VELP (IC-1. J) 

C01=C01* AZFAC 
CP1=CP1*REN0RM 
RCOS=RCOSP 
IX=IX-2 
GO  TO  134 


325  CONTINUE 
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32 


40 


FILE  2 

SO ( IC) =-l . + ( 1 . +S01 ) * ( 1 . +S02) *C02/ (C01* ( S02-S01 ) +C02* ( 1 . +S01 ) ) 

SP ( IC) =-l . + ( 1 . +SP1 ) * ( 1 . +SP2 ) *CP2/ (CPI* (SP2-SP1 ) +CP2* ( 1 . + SP1 ) ) 

DO  32  J=1 , 21 

VEL(IC, J) =VELO ( IC , J) +VELP ( IC , J) 

RVCT ( J) =VELP ( IC , J) 

IF ( IPRINT . NE . 4 ) GO  TO  1000 
WRITE (6.1) 

IC, SEN, ( CRRNTO ( IC , KX ) ,KX=1.5) 

IC  , RMTX (1,1) , AL , A , SO ( IC ) , SP ( IC ) 

IEL, (VELO(IC. J) , J=l,ll) 


WRITE (6, 3) 

WRITE (6, 3) 

WRITE ( 6 , 2) 

IND=1 

WRITE (6,2) 

IND  = 2 

WRITE (6,2) 

1000  RETURN 
END 

SUBROUTINE  BTODR ( DLNGTH 


IND.  (VELPdC.  J)  , J=l,ll) 


IND, (RVCT ( J) . J=l.ll) 


IC) 


COMMON  DIAMS ( 18 ) , DUCTO ( 3 , 17 , 6 ) , DUCT1 ( 3 , 17 . 6 ) , DUCT2 (3,17.6) .IEL.LP, 
X BEND ( 10 , 6 ) , SGRUPS ( 23 ) , DWTS ( 21 ) , CRRNTO ( 9 , 10 ) .RCRRNT ( 9 . 10 ) 

COMMON  RVCT (21) , RMTX (21. 21) .SO (9) ,SP(9) ,DZR0(9) ,NEXT(9) ,COF(9) , 

X ECC ( 9 ) , DELR ( 9 ) ,SNYNEU(23) ,SINVAL(10) ,C01.C02,VEL(9, 21) .IPRINT 
COMMON  RCOSO , RCOSP . RCOS . DGRUPS ( 23 ) . VELO (9.21), VELP (9,21) 

X , VEL1 ( 5 ) , SCALE . SWTS (23.21). WTSD (23.21) 

COMMON  IDOSE , HNDRSN ( 23 ) . SEPI ( 23 ) . STH ( 23 ) , SINR ( 23 ) , THMLN ( 23 ) 

X. ENGAM ( 21 ) . DIMEN, SCOBLT( 23) , AWTS(6,5) ,SGAM(23) , DGAM ( 23 ) ,AMOD(5. 5) 
DIMENSION  STRVEL(22) .TVLO (22) .TVLP (22) .SIDE (22, 22) .BACK (22. 22) . 

X VPB ( 22 ) . VPS ( 22 ) , VOS ( 22 ) , VOB ( 22 ) 

1 FORMAT (1H  ) 

2 FORMAT ( 13 , 11 ( 1PE10 .2)) 

3 FORMAT ( 13 , 9E12 . 4 ) 


STRSO=SO ( IC-1 ) 

DLBND=1 . 

DISTNS=DLNGTH+DLBND 
DO  10  1=1.21 

10  STRVEL ( I ) =VELO ( IC-1 , I ) 
TVO=0 . 

TVP=0 . 

C 

C DIRECT  COMPONENT 

CALL  DCTSLT ( DLNGTH. IC) 

DO  20  1=1,21 
TVLO (I ) =VELO ( IC , I ) 

TVLP ( I ) =VELP ( IC , I ) 
TVO=TVO+VELO (IC , I ) 

20  TVP=TVP+VELP ( IC , I ) 

RCO=TVO*2 . / ( 1 . 001+S0 ( IC) ) 
RCP=TVP*2 ./ (1.001+SP(IC) ) 


C 

C SIDE  COMPONENT 

DO  200  1=2,21 
200  VELO (IC-1, I )=0. 

VELO ( IC-1, 1)=1. 

VELP ( IC, 1 ) =0 . 

SO(IC-l ) =0 . 5 

CALL  DCTSLT ( DISTNS, IC) 

DO  205  1=1,20 

205  SIDE(1,I)=(VELP(IC,I+1)-.54*VELP(IC.I) )/.1162 
DO  210  1=2,21 
IP=I-1 

DO  207  J=1.IP 
207  SIDE ( I , J) =0 . 

DO  210  J=I , 21 

210  SIDE ( I , J) =SIDE (1-1 , J-l) 

VPS(1)=0. 

VOS(1)=0. 

DO  220  J=2 , 20 

VOS ( J) =SIDE ( J . J) * ( VELP ( IC-1 , J) - . 2055*STRVEL ( J-l ) ) 
VPS(J)=0. 

JP=J-1 

DO  220  1=1, JP 

VPS ( J) =VPS ( J) +SIDE (I , J) * VELP (IC-1. I) 

IF ( I . GT . 1 ) VPS( J)=VPS( J) -SIDE (I . J) * . 2055*STRVEL (1-1) 
220  CONTINUE 
SVO=0. 

SVP=0. 

DO  230  1=2.20 

TVLO ( I ) =TVLO (I ) + . 829*VOS ( I ) 

TVLP ( I ) =TVLP ( I ) + . 829*VPS ( I ) 

SVO=SVO+ . 829*V0S (I ) 

230  SVP=SVP+ . 829*VPS ( I ) 

TVO=TVO+SVO 

TVP=TVP+SVP 
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FILE  2 


RC0=RC0+SV0*2 . / ( 1 . 001+SP ( IC ) ) 
RCP=RCP+SVP*2 . / ( 1 . 001+SP ( IC) ) 
C 

C BACK  COMPONENT 

DO  300  1=2,21 
VELO ( IC-1 , I ) =0 . 

300  RVCT ( I ) =0 . 

RVCT ( 1 ) = 1 . 

VELO ( IC-1 , 1 ) =1 . 

CALL  DCTISO ( DISTNS , IC ) 

BVP=0 . 

BV0=0 . 

DO  305  1=1,21 
BVO=BVO+VELO ( IC . I ) 
BVP=BVP+VELP ( IC , I ) 

305  BACK (1,1) =VEL ( IC , I ) 

RC0B=2 . / ( 1 . 001+SO( IC) ) 

RCPB=2 ./( 1 . 001+SP ( IC) ) 

C 


SO(IC-l) =-0. 5 

CALL  DCTSLT ( DISTNS, IC) 

BVPM=0. 

DO  310  1=2,21 
BVPM=BVPM+VELP ( IC , I ) 

310  BACK ( 1 , 1 ) = ( BACK ( 1 „ I ) - 1 . 3 7 * VELP ( IC , I ) )/.205 
BACK (1,1)=BACK( 1,1)/. 205 

RCSPB=BVP*RCPB-1 . 4*BVPM*2 . / ( 1 . 001+SP ( IC) ) 
BVP=BVP-1.4*BVPM 
RCSPB=RCSPB/BVP 
C 

DO  315  1=2,21 
IP=I-1 

DO  312  J=1,IP 
312  BACK ( I , J) =0 . 

DO  315  J=I , 21 

315  BACK(I, J) =BACK(I-1, J-l) 

VOB ( 1 ) =0 . 

VPB(l) =0. 

DO  320  J=2 , 21 

VOB( J)=BACK( J, J)* (VELP (IC-1, J) +STRVEL ( J-l ) ) 

VPB ( J ) =0. 

JP=J-1 

DO  320  1=1, JP 

IF(I.GT.l)  VPB ( J) =VPB ( J) +BACK ( I , J) *STRVEL ( 1-1 ) 
320  VPB ( J) =VPB ( J) +BACK ( I , J) *VELP ( IC-1 , I ) 

SV0=0. 

SVP=0 . 

DO  330  1=1,21 
SVO=SVO+ .171*VOB(I) 

SVP=SVP+ . 171*VPB ( I ) 

VELO ( IC , I ) =TVLO ( I ) + . 171*VOB ( I ) 

VELP ( IC , I ) =TVLP ( I ) + , 171* VPB ( I ) 

330  VEL ( IC , I ) =VELO ( IC , I ) +VELP ( IC , I ) 

TVO=TVO+SVO 

RCO=RCO+SVO*RCOB 


TVP = TVP  + S VP 
RCP  = RCP  + SVP  *RCSPB 
SO(IC) =-l. +2, *TVO/RCO 
SP ( IC) =-l . +2 . *TVP/RCP 
VELP ( IC , 1 ) =0 . 

VEL ( IC , 1 ) =VELO ( IC , 1 ) 

DO  50  1=1,21 

VELO ( IC-1 , I ) =STRVEL ( I ) 

50  RVCT ( I ) =VEL ( IC , I ) 

SO ( IC-1 ) =STRSO 
IF ( IPRINT . NE . 4 ) GO  TO  1000 
41  WRITE (6,1) 

IC,S0(IC) ,SP(IC) , TVO , RCO , TVP , RCP , VELO ( IC , 2 ) ,VELP(IC,2 
I EL, (VELO(IC. J) , J=l,21,2) 

IEL. (VELP ( IC , J) , J=l,21,2) 


WRITE (6,3) 
WRITE (6, 2) 
WRITE (6, 2) 
IND=17 
WRITE (6, 2) 
WRITE (6, 2) 
1000  RETURN 
END 


IND. 

IND, 


(TVLO(J) , J=1.21,2) 
(TVLP(J) , J=1 , 21,2) 
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1:* 

2: 

3: 

4: 

5: 

6: 
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8: 
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SUBROUTINE  ROOMX ( IC , SVELO , SVELP ) 

COMMON  DIAMS ( 18 ) . DUCTO ( 3 , 17 , 6 ) , DUCT1 ( 3 , 17 . 6 ) , DUCT2 (3,17,6) . IEL.LP, 
X BEND (10.6), SGRUPS ( 23 ) . DWTS ( 21 ) , CRRNTO (9.10). RCRRNT (9,10) 

COMMON  RVCT ( 21 ) . RMTX ( 21 , 21 ) , SO ( 9 ) , SP ( 9 ) , DZRO ( 9 ) , NEXT ( 9 ) . COF ( 9 ) , 

X ECC ( 9 ) , DELR ( 9 ) , SNYNEU ( 23 ) , SINVAL ( 10 ) , C01 . C02 . VEL (9.21), IPRINT 
COMMON  RCOSO . RCOSP . RCOS . DGRUPS ( 23 ) . VELO (9.21), VELP (9,21) 

X , VEL1 ( 5 ) . SCALE , SWTS (23.21), WTSD (23,21) 

COMMON  XDOSE . HNDRSN ( 23 ) . SEPI ( 23 ) . STH ( 23 ) , SINR ( 23 ) , THMLN ( 23 ) 

X, ENGAM ( 21 ) .DIMEN. SCOBLT( 23) . AWTS(6,5) ,SGAM(23) ,DGAM(23) ,AMOD(5,5) 
DIMENSION  D ( 20 ) . DL ( 20 ) , W ( 20 ) , OVEL (22.20), PVEL (22,20), SVELO (22,20). 
* SVELP (22, 20) 

1 FORMAT (1H  ) 

2 FORMAT ( 13 , 11 ( 1PE10 . 2 ) ) 

3 FORMAT ( 11 ( 1PE10 ,2)) 

WTH=C0F ( IC) 

HT=ECC(IC) 

XL = DELR ( IC) 

AEX= SCALE 
AEXR=DZRO ( IC+1 ) 

IF (AEXR . GT . 1 . ) AEXR=1 . 

WD=DZRO ( IC) 

RMCASE=DZRO ( 1 ) 

WDP=COF ( IC+1 ) 

NW= (ECC ( IC+1 ) +1 . 01) 

NL=( DELR (IC+1) +1.01) 

NT=NW*NL 

DELW = WDP /ECC ( IC  + 1 ) 

DELL=XL/DELR ( IC+1 ) 

TN=SQRT ( 1 . -SO ( IC-1 ) **2) /SO(IC-l) 

A=SP ( IC-1 ) * (1 . +SO ( IC-1 ) ) / ( 1 . +SP ( IC-1 ) ) /SO ( IC-1 ) 
RR=AEX/(WTH*HT+HT*XL+XL*WTH) 

AL= . 75 

ALR=AL* ( 1 . -RR/2 . ) 

TNR=SQRT ( HT*WTH/3 . 14159 ) /XL 

POSITION  1,N=0  (UNREFLECTED,  IN  FRONT  OF  DUCT) 

OVEL (1,1)  =2 . / ( 1 . +SO ( IC-1 ) ) *AEXR 
PVEL ( 1 , 1 ) =2 . / ( 1 . +SP ( IC-1 ) ) *AEXR 

POSITIONS  2 THRU  NT.  N IS  2 OR  GREATER 
DO  11  1=2, NT 

OVEL ( 3 . I ) = ( 1 . +SO ( IC-1 ) ) *RR*ALR* *2 
PVEL ( 3 , I ) = ( 1 . +SP ( IC-1 ) ) *RR*ALR**2 
DO  11  J=4 . 21 

OVEL ( J. I ) =ALR*OVEL ( J-l , I ) 

11  PVEL ( J , I ) =ALR*PVEL ( J-l , I ) 

POSITION  1.  N IS  2 OR  GREATER 
DO  12  J=3 , 21 
OVEL ( J , 1 ) = . 5 *OVEL ( J , 2 ) 

12  PVEL ( J , 1 ) = . 5*PVEL ( J , 2 ) 

POSITIONS  2 THRU  NW,  N=0  (UNREFLECTED) 

DO  8 1=2, NW 
OVEL ( 1 , I ) =0 . 

8 PVEL ( 1 , I ) =0 . 

POSITIONS  ( NW+1 ) THRU  NT,N=0  (UNREFLECTED) 

DO  13  1=1, NL 
FI= ( 1-1 ) 

DO  13  J=1 , NW 
F J= ( J-l ) 

K= ( 1-1 ) *NW 

DL ( K+ J) =DELL*FI 

W ( K+ J) =DELW*FJ 

13  D ( K+ J) =SQRT ( DL (K+J)**2+W(K+J)**2  +.0001) 

K=NW+1 

DO  25  I = K , NT 

FR= ( DL ( I ) *TN-W ( I) + . 5*WD ) /WD 
IF(FR)  16.16,17 

16  FR=0. 

17  IF ( FR-1 . ) 19.18,18 

18  FR=1 . 

19  SAF=SQRT ( 1 , +AEX/3 .14159/D(I)**2) 

SAF=DL ( I ) /D ( I ) * ( 1 . -1 . /SAF) 

IF ( SAF- ( 1 . -SO ( IC-1 ) ) >20.20,21 
21  SAF=1 . -SO ( IC-1 ) 

20  OVEL (1,1) =FR*SAF/ ( 1 . -SO ( IC-1 ) ) 

PVEL ( 1 , I ) = ( 1 . +SP ( IC-1 ) ) * ( A*OVEL ( 1 , I ) / ( 1 . +SO ( IC-1 ) ) + ( 1 . -A) *SAF) 

25  CONTINUE 

MODIFY  FOR  BEND  INTO  ROOM  CASES 
IF (RMCASE)  100.51,150 

ENTERING  UNREFLECTED  MOVING  TOWARDS  DETECTORS 

100  K=NW+1 

DO  102  I=K , NT , NW 
OVEL ( 1 , I ) =0 . 

DO  102  J=2 , NW 


76 


FILE  3 


Continued 


83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 
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111 
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113 

114 
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123 

124 
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102  OVEL ( 1 , 1+ J-l ) =2 . *OVEL ( 1 , 1+ J-l ) 

GO  TO  51 

C ENTERING  UNREFLECTED  MOVING  TOWARDS  WALL 

150  DO  152  1=1. NT 
TEMP=DL ( I ) 

DL ( I ) =W ( I ) + ( WTH-WDP) 

W ( I ) =TEMP 

152  D ( I) =SQRT (DL( I )**2+W(I)**2+ .0001) 

TEMP=NW 

NW=NL 

NL=TEMP 

K=NW+1 

DO  154  I=K , NT 
154  OVEL ( 1 . 1 ) = 0 . 

GO  TO  32 

C ONCE-REFLECTED  COMPONENT  (N=l) 

51  SPX=1. 

RMX=0 . 

IF(TNR.GT.TN)  GO  TO  30 
SPX=0. 

RMX=1 . 

30  NLHF= (NL/2+1) 

DO  31  1=1 , NLHF 
KI= (1-1) *NW 
KIR= (NL-I ) *NW 
DO  31  J=1 , NW 
TEMP=DL( J+KI ) 

DL ( J+KI ) =DL ( J+KIR) 

31  DL( J+KIR) =TEMP 
DO  34  1=1. NT 

34  D ( 1 ) =SQRT(DL (I)**2+W(I)**2  +.0001) 

K=NT-NW+1 

DO  33  J=1 , NT 
IF(J.LT.K)  GO  TO  35 

OVEL (2. J) =OVEL ( 1 , J) *AL*2 . *DL ( J) /D( J) +RMX* ( 1 . +SO ( IC-1 ) ) *AL*RR/2 . 
PVEL ( 2, J) =PVEL ( 1 , J) *AL*2 . *DL ( J) /D ( J) +RMX* ( 1 . +SP ( IC-1 ) )*AL*RR/2. 
GO  TO  33 

35  SAFP=SQRT ( 1 . +AEX/3 .1415926/D(J)**2) 

3AFP=DL (J)/D(J)*(1.-1. /SAFP ) 

OVEL ( 2 , J) = ( 1 . +SO ( IC-1 ) ) *AL* ( SPX*SAFP+RMX*RR) 

PVEL ( 2 , J) = ( 1 . +SP ( IC-1 ) ) *AL* ( SPX* (A* SAFP/ (l.+SO(IC-l) )+() . -A) *RR) 

* +RMX*RR) 

33  CONTINUE 

GO  TO  57 

C RMCASE=1 . (N=l ) 

32  OVEL (2,1) =AL* ( 1 . +S0 ( IC-1 ) ) 

PVEL (2.1)  =AL* ( 1 . +SP ( IC-1 ) ) 

DO  36  J=2 , NT 

SAFP=SQRT (1.+AEX/3. 1415926/D (J)**2) 

SAFP=DL ( J) /D ( J) * ( 1 . -1 . /SAFP ) 

OVEL ( 2 , J) =AL* (l.+SO(IC-l) )*SAFP 

36  PVEL ( 2 , J) =AL* ( 1 . +SP ( IC-1 ) ) *SAFP 
TEMP=NW 

NW=NL 

NL=TEMP 

C RENORMALIZATION 

57  DO  55  1=1, NT 
DO  55  J=1 . 21 
IJ=I+J 

IF ( I J . EQ . 2)  GO  TO  55 

OVEL ( J, I ) =OVEL ( J. I ) *OVEL (1,1) 

PVEL ( J, I ) =PVEL ( J , I ) *PVEL (1,1) 

55  CONTINUE 

C MATRIX  GENERATION  AND  MULTIPLICATION ( S ) 

DO  70  1=1. NT 
DO  58  J=1 , 21 
SVELP ( J , I ) =0 . 

58  SVELO ( J, I) =OVEL (1,1) *VELO ( IC-1 . J) +PVEL (1,1) *VELP ( IC-1 . J) 

DO  65  IP=2 , 21 

DO  65  J=2 . IP 

65  SVELP ( IP, I ) = SVELP ( IP , I ) +OVEL ( J, I) *VELO ( IC-1 , IP+1- J) +PVEL ( J, I ) * 

* VELP (IC-l.IP+l-J) 

70  CONTINUE 

IF ( IPRINT . NE . 4 ) GO  TO  1000 
40  WRITE (6,1) 

WRITE (6,2)  IC . A EX , AEXR . WD , RMCASE . DELW , DELL , TN , TNR , A , RR , ALR 
WRITE (6, 1) 

WRITE (6,3)  ( SVELO (J,l) . J=1 , 21 , 2) 

WRITE (6. 3)  ( SVELO ( J, 2 ) , J=1 ,21,2) 

WRITE (6,3)  ( SVELO ( J , NT ) , J=1 , 21 , 2 ) 

WRITE (6,3)  (OVEL(J.l) . J=l,10) 

WRITE (6, 1) 

WRITE (6, 3)  ( SVELP (J.l) . J=l,21,2) 
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( SVELP ( J , 2 ) , J=l,21.2) 
(SVELP(J.NT) . J=1.21.2) 
( PVEL (J,l) , J=1 , 10 ) 


WRITE (6. 3) 

WRITE (6, 3) 

WRITE (6. 3) 

WRITE (6.1) 

1000  RETURN 
END 

SUBROUTINE  DOSX ( IC . SVELO , SVELP ) 

COMMON  DIAMS ( 18 ) . DUCTO ( 3 , 17 , 6 ) . DUCT1 ( 3 . 17 . 6 ) , DUCT2 (3.17,6) .IEL.LP. 
X BEND ( 10 , 6 ) , SGRUPS ( 23 ) , DWTS ( 21 ) , CRRNTO (9,10) , RCRRNT (9,10) 

COMMON  RVCT ( 21 ) , RMTX ( 21 . 21 ) , SO ( 9 ) , SP ( 9 ) , DZRO ( 9 ) , NEXT ( 9 ) . COF ( 9 ) , 

X ECC ( 9 ) , DELR ( 9 ) , SNYNEU ( 23 ) , SINVAL ( 10 ) . COl . C02 . VEL (9,21), IPRINT 
COMMON  RCOSO . RCOSP . RCOS , DGRUPS ( 23 ) , VELO (9.21). VELP (9,21) 

X , VEL1 ( 5 ) , SCALE , SWTS (23,21) , WTSD (23,21) 

COMMON  IDOSE . HNDRSN ( 23 ) , SEPI ( 23 ) , STH ( 23 ) . SINR ( 23 ) , THMLN ( 23 ) 

X , ENGAM ( 21 ) , DIMEN. SCOBLT( 23) . AWTS(6.5) ,SGAM(23) ,DGAM(23) ,AMOD(5.5) 
DIMENSION  SVELO (22. 20) , SVELP ( 22 . 20 ) .PRAT (23. 2) ,ORAT(20,2) . 

* DRAT (20.2). SPCTRM ( 23 ) . DOSE ( 20 ) , DFLXN ( 21 ) 

1 FORMAT (1H  ) 

2 FORMAT (5H  IEL=.I2.5H  IC=.I2) 

3 FORMAT (8H  DIR/ICD) 

4 FORMAT (8H  REF/ICD) 

5 FORMAT (8H  TOT/ICD) 

6 FORMAT (8H  DIR/IFD) 

7 FORMAT (8H  REF/IFD) 

8 FORMAT (30H  DOSE/ENTRANCE  DOSE 

10  FORMAT ( 9 ( 1PE12 . 4 ) ) 

11  FORMAT (8H  RC0SIN= , 1PE10 . 3 , 7H 
X7H  EXF= , 1PE10 . 3 ) 

12  FORMAT ( 55H)  MEG, N 
FORMAT (5H  DOSE) 

FORMAT ( 55H  VEL(N)/IC=1 
FORMAT ( 15 . 7 ( 1PE10 . 3 ) ) 

FORMAT (5H  SO=  ,7(1PE10 
FORMAT (5H  SP  = . 7 ( 1PE10 
RCOSIN=2 . / ( 1 . + COF ( 1 ) ) 

CDIN=DWTS ( 1 ) 

IF (CDIN . EQ . 0 . ) CDIN=1 . 

NW= ( ECC (IC) +1.01) 

NL=  ( DELR  (10+1.01) 

NT=NW*NL 
NMAX--  21 

IF ( IDOSE. NE. 5)  GO  TO  122 
NMAX=10 

122  DO  100  MEG=1 . 23 
SPCTRM (MEG) =0 . 

U=l. 

DO  100  N=1 . NMAX 

SPCTRM ( MEG ) = SPCTRM ( MEG ) + ( SVELO (N, 1) + SVELP (N.1))*SWTS(MEG.N)/U 
100  U=U* . 696 


W ( FT) = , F7 . 2 , 9H  L(FT)=.F7.2) 


CDIN= , 1PE10 . 3 , 7H  DIN= , 1PE10 . 3 . 
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SOURCE 


3)  ) 
.3)  ) 


SPCTRM 


DFLXN 


VELO 


4 


VELP  ) 


DO  50  1=1. NT 
ORAT ( 1 , 1 ) = 0 . 

PRAT ( I , 1 ) =0 . 

U=i. 

DO  40  NREF=1 , NMAX 

ORAT (1,1) =ORAT (1,1) + SVELO ( NREF . I ) *DWTS ( NREF ) /U 
PRAT ( 1 , 1 ) = PRAT ( 1 , 1 ) + SVELP ( NREF . I ) * DWTS ( NREF ) /U 
DFLXN ( NREF ) = ( SVELP ( NREF . 1 ) + SVELO ( NREF , 1 ) ) *DWTS ( NREF ) /U 
40  U=U* . 696 

DRAT (1,1) =PRAT (1,1) +ORAT ( I . 1 ) 

DOSE ( I ) =DRAT (1,1) 

50  CONTINUE 

IF (NEXT ( 1 ) ) 101.101.102 
102  ZD=DELR ( 1 ) 

EXF= ( . 671 5+ ZD* ( . 924+ . 565*ZD) ) / ( . 0685+ZD* ( . 5445+ . 137*ZD) ) 
CDIN=CDIN*EXP ( ZD*EXF ) 

101  DIN=CDIN*RCOSIN 


DO  120  1=1, NT 
ORAT (1,2) =0RAT (1,1) /DIN 
PRAT (1.2) =PRAT (1,1) /DIN 
DRAT ( I , 2 ) = DRAT (1,1) /DIN 
ORAT (1,1) =ORAT (1,1) /CDIN 
PRAT (1,1) =PRAT (1,1) /CDIN 
120  DRAT (1,1) =DRAT (1,1) /CDIN 
IF ( IPRINT. EQ.O)  GO  TO  1000 
IF ( IPRINT . LE . 2 ) GO  TO  128 
WRITE (6. 2)  I EL , IC 
105  WRITE (6,1) 

WRITE (6.12) 

JC=IC-2 
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FILE  3 Continued 

DO  109  1=1.21 

WRITE (6,15)1, SGRUPS ( I ) .SPCTRM(I) .DFLXN(I) ,VELO( JC.I) ,VELP( JC,I) 

109  CONTINUE 
MEG=22 

WRITE (6.15)  MEG . SGRUPS ( MEG ) . SPCTRM ( MEG ) 

MEG=23 

WRITE (6,15)  MEG . SGRUPS ( MEG ) . SPCTRM ( MEG ) 

WRITE (6.1) 

WRITE (6.1) 

WRITE (6. 14) 

DO  113  N=1 , 21 

WRITE (6,15)  N. (VEL(ID.N) ,ID=1, JC) 

113  CONTINUE 

WRITE (6,16)  (SO(ID) ,ID=1, JC) 

WRITE (6,17)  (SP(ID) ,ID=1. JC) 

WRITE (6,1) 

IF (IPRINT.LT. 4)  GO  TO  128 
WRITE (6. 11)  RCOSIN,  CDIN.DIN.EXF 
WRITE (6, 3) 

WRITE (6,10)  ( ORAT (1,1) ,1=1, NT ) 

WRITE (6,4) 

WRITE (6. 10)  ( PRAT (1,1) .1=1, NT ) 

WRITE (6,1) 

WRITE (6. 5) 

DO  125  J = 1 , NW 

WRITE(6„ 10)  ( DRAT ( 1+ J-l , 1) ,I=1,NT,NW) 

125  CONTINUE 
WRITE (6. 6) 

WRITE (6. 10)  ( ORAT (1.2) ,1=1. NT) 

WRITE (6, 7) 

WRITE (6,10)  ( PRAT (1,2) ,1=1, NT ) 

WRITE (6,1) 

128  WRITE (6. 8)  COF ( IC ) , DELR ( IC-1 ) 

WRITE (6,1) 

DO  130  J=1 , NW 

WRITE (6. 10)  ( DRAT (1+ J-l. 2) ,I=1,NT.NW) 

130  CONTINUE 
WRITE (6.1) 

WRITE (6. 13) 

DO  135  J=1 , NW 

WRITE (6. 10)  ( DOSE (1+ J-l) ,I=1.NT,NW) 

135  CONTINUE 
1000  RETURN 
END 

L'O 

SUBROUTINE  DOSE(IC) 

COMMON  DIAMS ( 18 ) . DUCTO ( 3 , 17 , 6 ) , DUCT1 ( 3 , 17 , 6 ) , DUCT2 (3,17,6) ,IEL,LP, 

X BEND (10,6), SGRUPS ( 23 ) , DWTS ( 21 ) , CRRNTO (9,10). RCRRNT (9,10) 

COMMON  RVCT ( 21 ) , RMTX ( 21 , 21 ) , SO ( 9 ) , SP ( 9 ) , DZRO ( 9 ) . NEXT ( 9 ) . COF ( 9 ) , 

X ECC ( 9 ) , DELR ( 9 ) , SNYNEU ( 23 ) , SINVAL ( 10 ) , COl . C02 , VEL (9,21), IPRINT 
COMMON  RCOSO . RCOSP . RCOS . DGRUPS ( 23 ) , VELO (9,21), VELP (9.21) 

X . VEL1 ( 5 ) , SCALE . SWTS (23,21), WTSD (23,21) 

COMMON  IDOSE , HNDRSN ( 23 ) , SEPI ( 23 ) , STH ( 23 ) , SINR ( 23 ) , THMLN ( 23 ) 

X , ENGAM ( 21 ) , DIMEN . SCOBLT ( 23 ) , AWTS (6.5), SGAM ( 23 ) , DGAM ( 23 ) , AMOD (5.5) 

DIMENSION  VELX ( 21 ) , SPCTRM ( 23 ) , DFLXN ( 21 ) , DELO ( 9 ) . DELP ( 9 ) , 

X DOUT ( 9 ) , REDUXT ( 9 ) , REDUXD ( 9 ) 

1 FORMAT (1H  ) 

2 FORMAT (1H1) 

3 FORMAT ( 15 , 7 ( 1PE11 . 3 ) ) 

4 FORMAT (55H  IC  DELP  DELO  DELX  RCOSO  RCOSP  RCOSX  ) 

5 FORMAT ( 5 5H  NEXT(JC)  CDIN  CDOUT  DIN  DOUT  REDUXT (JC)  REDUXC  ) 

7 FORMAT (9H  NEXT ( I ) = , 12 . 7H  DOSE= , 1PE10 . 3 . 14H  REDUX  FACTR= , 1PE10 . 3 ) 

8 FORMAT ( 15 , 5 ( 1PE10 . 3 ) ) 

9 FORMAT ( 5 5H  MEG, N SOURCE  SPCTRM  DFLXN  VELO  VELP  ) 

10  FORMAT ( 15 , 4 ( 1PE10 . 3 ) . 110 . 4 ( 1PE10 . 3) ) 

11  FORMAT ( 55H  VEL(N)/IC=1  2 3 4 5 ) 

12  FORMAT ( 15 , 10 ( 1PE10 . 3 ) ) 

13  FORMAT (5H  SO=  . 10 ( 1PE10 . 3) ) 

14  FORMAT (5H  SP=  , 10 ( 1PE10 . 3 ) ) 

15  FORMAT (7H  ARTO=  .1PE10.3) 

RCOSO=2 . / ( 1 . +SO ( IC-1 ) ) 

RCOSP=2 . / ( 1 . +SP ( IC-1 ) ) 

RCOSIN=2 . / ( 1 . +COF ( 1 ) ) 

CALL  RETRN ( VELX , SX , IC) 

RCOSX=2 . / (1 . +SX) 

CDIN=0 . 

JC=IC-1 
DO  19  1=1, JC 
DELO ( I ) =0 . 

19  DELP ( I ) =0 . 

DELX=0 . 

CDR1=0 . 

NMAX=21 
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331 

IF(IDOSE. LT. 5)  GO  TO  20 

332 

NMAX=10 

333 

20 

CDIN=DWTS ( 1 ) 

334 

IF (CDIN . EQ . 0 . ) CDIN=1 . 

335 

U=l. 

336 

CRD1=CRD1+VEL ( IC-1 , 2 ) *DWTS ( 2 ) /U 

337 

DO  100  NREF = 1 . NMAX 

338 

DO  99  1 = 1.  JC 

339 

DELO ( I ) =DELO ( I ) +VELO ( I . NREF ) *DWTS ( NREF ) /U 

340 

99 

DELP ( I ) =DELP ( I ) +VELP ( I , NREF ) *DWTS ( NREF ) /U 

341 

DELX=DELX+ VELX ( NREF ) * DWTS ( NREF ) /U 

342 

100 

U=U*0.6y6 

343 

DELX=DELX* SCALE 

344 

DO  213  M=1 . 23 

345 

SPCTRM(M) =0 . 

346 

U=l. 

347 

V=l. 

348 

DO  212  N=1 , NMAX 

349 

DFLXN (N) = ( VELO ( JC . N ) *RCOSO+VELX ( N ) *RC0SX ) *V/U 

350 

SPCTRM (M) =SPCTRM (M) + ( SWTS (M.N)*1.E10) *DFLXN (N ) /V 

351 

V=V* . 75 

352 

212 

U=U* . 696 

353 

213 

SPCTRM (M) =SPCTRM(M) /I .E10 

3 54 

IF (NEXT (1) ) 101.101,102 

355 

C 

3 56 

102 

ZD=DELR ( 1 ) 

357 

EXF= ( ,6715+ZD*( .924+.565*ZD) )/( ,0685+ZD*( 

. 5445+ . 137*ZD ) ) 

358 

CDIN=CDIN*EXP ( ZD*EXF ) 

359 

C 

360 

101 

DIN=CDIN*RCOSIN 

361 

CDOUT=DELO ( JC ) +DELX 

362 

DO  103  1=1, JC 

363 

DOUT  ( I ) =DELO  (I)  *2  . / (1 . +SO  (I)  ) -s-DELP  (I)  *2 . / (1 . +SP(I)  ) 

3M 

REDUXT ( I ) =DOUT ( I ) /CDIN 

3^5 

103 

REDUXD ( I ) =DOUT ( I ) /DIN 

3t  6 

REDUXC=CDOUT/CDIN 

3 1 7 

IF ( IPRINT . EQ . 0 ) GO  TO  1000 

36  8 

IF ( IPRINT . EQ . 1 ) GO  TO  106 

36  9 

DO  104  1=1, JC 

370 

WRITE (6,7)  NEXT ( I ) , DOUT ( I ) . REDUXD ( I ) 

371 

104 

CONTINUE 

372 

GO  TO  107 

373 

106 

WRITE (6, 7)  NEXT( JC) ,DOUT( JC) , REDUXD (JC) 

374 

107 

IF (NEXT ( JC+2 ) . EQ.8)  GO  TO  1000 

375 

WRITE (6, 1) 

376 

IF ( IPRINT. GT. 2)  GO  TO  105 

3.  7 

GO  TO  1000 

378 

105 

WRITE (6,1) 

3/9 

WRITE (6, 15)  SCALE 

360 

WRITE (6. 9) 

381 

DO  109  1=1.21 

382 

WRITE (6,8)1, SGRUPS ( I ) , SPCTRM ( I ) . DFLXN ( I ) , 

VELO ( JC , I ) , VELP ( JC , I ) 

383 

109 

CONTINUE 

384 

MEG=22 

385 

WRITE (6.8)  MEG , SGRUPS ( MEG ) . SPCTRM ( MEG ) 

386 

MEG=23 

387 

WRITE (6,8)  MEG . SGRUPS ( MEG ) . SPCTRM ( MEG ) 

388 

WRITE (6.1) 

389 

WRITE (6. 11) 

390 

DO  113  N=1 , 21 

391 

WRITE (6,12)  N. (VEL(ID.N) ,ID=1, JC) 

392 

113 

CONTINUE 

393 

WRITE (6,13)  ( SO ( ID) , ID=1 , JC) 

394 

WRITE (6,14)  (SP(ID) ,ID=1. JC) 

395 

110 

WRITE (6,1) 

396 

IF ( IPRINT. LE. 3) GO  TO  1000 

397 

WRITE (6, 4) 

398 

WRITE (6,3)  IC , DELP . DELO . DELX , RCOSO . RCOSP , 

RCOSX 

399 

WRITE (6. 5) 

400 

WRITE (6.3)  NEXT ( JC ) , CDIN , CDOUT . DIN , DOUT ( JC ) , REDUXT ( JC ) . REDUXC 

401 

1000 

RETURN 

402 

END 
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8: 
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SUBROUTINE  RETRN (VELX, SX, IC) 

COMMON  DIAMS ( 18 ) . DUCTO ( 3 . 17 , 6 ) . DUCT1 ( 3 . 17 . 6 ) . DUCT2 (3.17,6) .IEL.LP. 
X BEND (10.6). SGRUPS ( 23 ) . DWTS ( 21 ) . CRRNTO (9.10). RCRRNT (9.10) 
COMMON  RVCT ( 21 ) . RMTX ( 21 . 21 ) . SO ( 9 ) . SP ( 9 ) . DZRO ( 9 ) . NEXT ( 9 ) . COF ( 9 ) , 

X ECC ( 9 ) . DELR ( 9 ) . SNYNEU ( 23 ) , SINVAL ( 10 ) . COl . C02 . VEL (9.21). IPRINT 
COMMON  RCOSO . RCOSP . RCOS . DGRUPS ( 23 ) . VELO (9,21), VELP (9,21) 

X , VEL1 ( 5 ) . SCALE . SWTS (23,21), WTSD (23.21) 

COMMON  IDOSE , HNDRSN ( 23 ) , SEPI ( 23 ) , STH ( 23 ) . SINR ( 23 ) . THMLN ( 23 ) 

X . ENGAM ( 21 ) , DIMEN , SCOBLT ( 23 ) . AWTS (6,5), SGAM ( 23 ) . DGAM ( 23 ) . AMOD (5,5) 
DIMENSION  VELX ( 21 ) , VELPX ( 21 ) . SVELP ( 21 ) , SVELO ( 21 ) . XVELO ( 21 ) . 

X XVELP ( 21 ) . SVEL ( 21 ) 

1 FORMAT (1H  ) 

2 FORMAT ( 13 , 11 ( 1PE10 . 2 ) ) 

3 FORMAT ( 13 , 7E13 . 3 ) 

KSUP=DZRO ( IC ) 

DO  8 IR=1 , 21 

SVEL ( IR ) =VEL ( IC- 1 , IR ) 

VELPX ( IR ) = VELP ( IC-1 , IR) 

8 VELX ( IR )= VELP (IC-l.IR) 

SX=SP(IC-1) 

SXP=SP(IC-1) 

IF(KSUP)  50.50.10 
10  RDLNTH=COF ( IC) 

IF ( RDLNTH . GT . 5 . ) GO  TO  50 
DLNGTH=COF ( IC-1 ) 

IF ( RDLNTH. LT. DLNGTH + 1. ) GO  TO  50 
IF(RDLNTH.GT.DLNGTH+5. ) GO  TO  50 
TOTO=0 . 

T0TP=0 . 

DO  12  NREF=1 , 21 
TOTO=TOTO+VELO ( IC-1 . NREF) 

12  TOTP=TOTP+VELP( IC-1, NREF) 

RC= (TOTO*RCOSO+TOTP*RCOSP) / (TOTO+TOTP) 

IF (RC-1 . ) 13. 13. 14 

13  RC=1 . 01 

14  STOTO=0 . 

ST0TP=0 . 

XTOTO=0 . 

XT0TP=0 . 

SSP=0. 

XSP=0. 

XS0=0. 

SO(IC+1)=0. 

SP(IC+1)=0. 

26  DO  27  IR=1 , 21 

RVCT ( IR ) =VEL ( IC-2 , IR) 

SVELO ( IR )= VELO ( IC- 1 . IR ) 

27  SVELP (IR)=VELP( IC-l.IR) 

SSO=SO(IC-l) 

SSP=SP(IC-1) 

IEL=NEXT ( IC-1 ) 

LP  = 2 

GO  TO  (100.110,300.500,500) ,IEL 
100  CALL  DCTISO (RDLNTH, IC-1) 

GO  TO  30 

110  CALL  DCTSLT (RDLNTH, IC-1) 

GO  TO  30 
300  DLNGTH =0 . 

500  IF  ( IEL . GE . 4 ) CALL  BT0DA( RDLNTH . IC-1 ) 

30  RISO=RDLNTH-DLNGTH 
RSLT=RISO- . 75 
LP=3 

RVCT ( 1 ) =0 . 

DO  35  IR=1 , 20 

35  RVCT (IR+1)=0.53*VEL( IC-l.IR) 

CALL  DCTISO (RISO.IC) 

DO  40  IR=1 , 21 

XVELO  ( IR ) = VELO  ( IC . IR ) 

XVELP ( IR) = VELP (IC.IR) 

IF ( IR . GT . 20 ) GO  TO  41 

40  VELO  (IC-1,  IR+1)  = .06*VEL(  IC-1.  IR) 

41  XSO=SO(IC) 

XSP=SP(IC) 

SO(IC-1)=0. 

CALL  DCTSLT ( RSLT. IC) 

DO  45  IR=1 . 21 

VELO (IC-l.IR) = SVELO ( IR ) 

VELP ( IC- 1 . IR ) = SVELP ( IR ) 
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FILE  4 Continued 


83:  VEL(IC-1.IR)=SVEL(IR) 

84  : VELX ( IR ) =XVELO ( IR ) +XVELP ( IR ) + VELO ( IC , IR ) + VELP ( IC . IR ) +VELPX ( IR ) 

85:  XTOTO=XTOTO+XVELO ( IR) 

86:  XTOTP=XTOTP+XVELP(IR) 

87:.  STOTO=STOTO+VELO ( IC, IR) 

88:  45  STOTP=STOTP+VELP ( IC . IR) 

89:  SO(IC-l)=SSO 

90:  SP(IC-1)=SSP 

91 : TOTX= TOTP  + STOTO  + STOTP +XTOTO + XTOTP 

92  : RCOS  = 2 . * ( TOTP/ ( 1 . +SXP ) +XTOTO/ ( 1 . +XSO ) + XTOTP/ ( 1 . + XSP ) 

93:  X +STOTO/ ( 1 . +S0 ( IC) ) +STOTP/ ( 1 . +SP ( IC) ) ) /TOTX 

94:  SX=2 . /RCOS-1 . 

95:  C 

96:  IF ( IPRIMT . NE . 4 ) GO  TO  50 

97:  47  WRITE (6,1) 

98:  WRITE (6.2)  IEL. (XVELO(IR) . IR=1« 21. 2) 

99:  WRITE (6,2)  IC, ( VELPX(IR) . IR=1 - 21 . 2) 

100:  WRITE (6,2)  IC, (XVELP (IR) , IR=1 , 21 . 2) 

101:  WRITE (6,2)  IC, (VELX(IR) . IR=1, 21, 2) 

102:  WRITE (6,3)  IC.  SX,  SSP, XSO, XSP , SO ( IC) , SP ( IC) , DLNGTH 

103:  WRITE (6, 3)  IC,  TOTP , XTOTO. XTOTP , STOTO , STOTP, TOTX .RDLNTH 

104:  WRITE (6.1) 

105:  50  RETURN 

106:  END 

107:  SUBROUTINE  WTLISTS  ( JSTT, JXWT. IOX, NMAX, DWTS , SWTD , ENGAM, IDOSE, 

108:  X ISORS.LP) 

109:  DIMENSION  DWTS (21) , ENGAM ( 21 ) , SWTD (21) . IXFILE ( 65 , 3 ) . WTFILE ( 65 . 21 , 3 ) 

110:  4 FORMAT (315. 3112) 

111:  5 FORMAT (31 5) 

112:  8 FORMAT ( 5 ( 1PE12 . 4 ) ) 

113:  C - - - - -SETS  UP  AND  READS  PRECALCULATED  WEIGHTS  FROM  A FILE  (WTSFILE) 

114:  ITOT=65 

115:  15  IF (IOX)  20.17.100 

116:  17  IF ( IXALT)  80,1000.50 

117:  C (MAIN  PROGT^AM:  IOX=ql) 

118:  C - - - - - SET  UP  3 INDEX  FILES  WITH  TIE  INTERNAL  WTFILES 

119:  20  DO  40  I=l.ITOT 

120:  READ (7.5)  ( IXFILE ( I . J)  . J=1 , 3) 

121:  DO  40  IX=1 , 21 

122:  READ (7,8)  (WTFILE ( I. IX, J) . J=1 . 3 ) 

123:  40  CONTINUE 

124:  IOX=Q 

125:  IXALT=1 

126 : REWIND  7 

127:  C - - WRITE (4, 4)  ( IXFILE ( 12 , J) , J=1.3) 

128:  C - - WRITE (4,8)  (WTFILE (12 , 2 , J) , J=1.3) 

129:  GO  TO  1000 

130:  C - - - - - WRITE  BACK  THE  INFORMATION  INTO  PREWTS 

131:  100  DO  110  1=1 , ITOT 

132:  WRITE (7.5)  ( IXF ILE(I.J) ,J=1,3) 

133:  DO  110  NX=1 , 21 

134:  WRITE (7.8)  ( WTFILE (I, NX. J)  . J=l,3) 

135:  110  CONTINUE 

136:  REWIND  7 

137:  GO  TO  1000 

138:  C (IOX  VALUE  OUTPUT  DOESN'T  MATTER) 

139:  C - - - - - NEXT  EITHER  RIAD  THE  OLD  WEIGHTS  OR  STORE  NEW  ONES 
140:  C ( JSTT=ql . CALC  THEN  WRITE:  JSTT=0=ISTT,  READ  THEN  SKIP) 

141:  50  IXALT=- IXALT 

142:  JXWT=IDOSE»5*ISORS+20*LP-25 

143:  IF ( JSTT)  1000,51.1000 

144:  51  JSTT= IXFILE ( JXWT, 2) 

145:  IF ( JSTT . NE . 0 ) GO  TO  1000 

146:  52  JMAX=IXFILE( JXWT. 3) 

147:  DO  55  NX=1.JMAX 

148:  DWTS ( NX ) =WTFILE ( JXWT . NX . 1 ) 

149:  SWTD ( NX ) =WTFILE ( JXWT . NX , 2 ) 

150:  55  ENGAM ( NX )= WTFILE (JXWT. NX. 3) 

151:  NMAX= JMAX 

152:  GO  TO  1000 

153:  C - - - - - WRITE  WEIGHTS  JUST  COMPUTED 
154:  80  IXALT  = - IXALT 

155:  IF  ( JSTT  . EQ  . 0 ) GO  TO  1000 

156:  IXFILE ( JXWT. 2) =0 

157:  IXFILE ( JXWT. 3) =NMAX 

158:  DO  85  NX=1.NMAX 

159:  WTFILE ( JXWT . NX , 1 ) = DWTS ( NX ) 

160:  WTFILE (JXWT, NX, 2) = SWTD (NX) 

161:  WTFILE ( JXWT , NX , 3) = ENGAM (NX) 

162:  85  CONTINUE 

163:  1000  RETURN 

164:  END 

165:  SUBROUTINE  ENTDAT 
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COMMON  DIAMS ( 18 ) . DUCTO ( 3 . 17 . 6 ) , DUCT1 ( 3 . 17 , 6 ) . DUCT2 (3.17,6) .IEL.LP. 
X BEND (10.6). SGRUPS ( 23 ) , DWTS ( 21 ) . CRRNTO (9.10). RCRRNT (9.10) 

COMMON  RVCT ( 21 ) . RMTX ( 21 . 21 ) . SO ( 9 ) . SP ( 9 ) . DZRO ( 9 ) . NEXT ( 9 ) . COF ( 9 ) . 

X ECC ( 9 ) . DELR ( 9 ) , SNYNEU( 23) , SINVAL ( 10 ) , COl , C02 , VEL (9,21) . IPRINT 
COMMON  RCOSO . RCOSP , RCOS . DGRUPS ( 23 ) , VELO ( 9 . 21 ) . VELP ( 9 . 21 ) 

X , VEL1 ( 5 ) .SCALE. SWTS (23. 21) ,WTSD(23,21) 

COMMON  IDOSE , HNDRSN ( 23 ) . SEPI ( 23 ) , STH ( 23 ) . SINR ( 23 ) . THMLN ( 23 ) 
X,ENGAM(21) , DIMEN. SCOBLT ( 23) ,AWTS(6. 5) .SGAM(23) ,DGAM(23) . AMOD ( 5 . 5 ) 
DIMENSION  TPX1 ( 17) . TPX2(17).  TPX3(51),  TPX4(51),  TPX5(17). 

X TPX6 ( 51 ) , TPX7 ( 17 ) , TPX8(51).  TPX9(51),  TP10(17),  TP11(51). 

X TP12(17) . TP13 ( 51 ) , TP14 ( 51 ) , TP15(17).  TP16(51).  TP17(  8). 

X TP18 ( 8).  TP19 ( 8).  TP20 ( 8),  TP21 ( 8).  TP22(  8).  TP23(23). 

X TP24(23)  TP25 ( 23 ) , TP26(23).  TP27(20).  TP28(10).  TP29(23). 

X TP30 ( 23 ) , TP31 ( 23 ) . TP32(23).  TP33(23) 

C 

DATA  TPX1  /0... 5,1. ,1.5. 2. .2. 5. 3. .3. 5. 4. .4. 5. 5. .5. 5. 

X 6. .6.5.7. .7.5.8./ 

DO  101  IL=1 , 17 
101  DIAMS (IL)=TPX1(IL) 

C 

DATA  TPX2  /l. , .413, .205. .111. .0685. .0481. .0354. 

X .0292, .0226, .0191. .0155. .0127. .0108. .0091. .0077. .0065. .0055/ 

DATA  TPX3  /. 35, .56. .66, .71, .74. .75. .76. 

X .753. .745. .73. .715, .698. .675. .655. . 635 . . 615 , . 60 . . 35 . .488. .595. 

X .675. .72. .747. . 76 . . 753 . . 745 . . 73 , . 715 , .698. . 675 . . 655 , . 635 . .615. 

X .60, .52, .555. .61, .725. .77. .79. .804. .80. .79. .775. .758. .733. .71. 

X .685. .66. .64. .62/ 

DATA  TPX4  /. 30. .59. .75. .85. .93, .98.1.02. 

X 1.04.1.05.1.05.1.05,1.05.1.05,1.05,1.05,1.05.1.05. .30. .705. .835, 

X .91, .975.1.02,1.05.1.06.1.07,1.08,1.08.1.08,1.08,1.08,1.08.1.08. 

X 1.08. .74. .825. .90. .965.1.015.1.05,1.06,1.07.1.08,1.08.1.08.1.08. 

X 1.08.1.08.1.08,1.08.1.08/ 

DATA  TPX5  /l. . .517, .285. .154. .097. .063. .046. .037. 

X .028. .022. .019. .016. .013. .0108, .009. .0074. .0062/ 

DATA  TPX3  /. 265. .38, .48, .55, .605, .65. .69. 

X .725, .75. .78, .80, .823, . 84 , . 85 . . 86 . . 87 . . 88 . .265. . 365 . . 46 . . 53 , . 58 , 

X .625. .66  .693. .72. .745, . 765 . . 785 . . 8 . . 82 , . 83 , . 84 . . 85 . . 40 . . 40 . . 43 . 

X .49, .54,  584, .62, .65. .675, .70. .72. .74. .75. .768. .78. .79. .80/ 

IDX=0 

DO  106  J=1 , 3 
DO  106  IL=1 . 17 
IDX=IDX+1 

IF(J-l)  106.104.105 

104  DUCTO ( 1 . IL , 1 ) =TPX2 ( IDX ) 

DUCTO ( 1 , IL . 4 ) =TPX5 ( IDX) 

105  DUCTO (J, IL, 2) =TPX3 (IDX) 

DUCTO ( J . IL , 3 ) =TPX4 ( IDX ) 

DUCTO ( J . IL . 5 ) =TPX6 ( IDX ) 

106  CONTINUE 
C 

DATA  TPX7  /2..1.,1.,1.,1.,1.,1.,1..1.,1.,1.,1., 

X 1..1..1..1..1./ 

DATA  TPX3  /. 35. .62. .52, .51. .51, .512, .518, 

X .523. .533. .548, .562, .58, .6. .622, .645, . 665 , . 685 . . 35 . . 59 . . 515 , . 51 . 

X .51, .512. .518, .523, .533. .548, .562. .58, .6. .622, . 645 . . 665 , . 685 , 

X .615. .56. .52. .51. .51. .512, .518. .523, .533. .548. .562. .58. .6. .622. 

X .645, .665. .685/ 

DATA  TPX9  / . 3 . . 61 . . 83 . . 945 , 1 . 00 , 1 . 04 . 1 . 06 
X. 1.065, 1.067. 1.064. 1.06. 1.04. 1.02. 1.00.  . 98 . . 96 . . 94 . . 3 . . 68 . . 865 . 

X .965.1.025.1.06.1.08.1.084.1.08.1.07,1.05.1.04.1.02.1.00. .98. .96. 
X .94, .68, .79, .905,1. .1.07.1.11,1.125,1.13.1.12.1.11.1.1.1.07.1.05. 
X 1.03.1.01. .99. .97/ 

DATA  TP10  /. 5.  .5. .5, .5.  .5. .5. .5. .5.  .5. .5, .5,  .5, 

X .5,  .5.  .5,  .5, .5/ 

DATA  TP11  /. 35.  .48. .58.  .64. .69, .735. .766, 

X .79. .81, .83, .843. .85. .86, .866. .872, .878, .884. .33.-44, .533. .59. .64 
X , .68, .71. .733. .75. .765. .775, .78. .786. .791. . 795 . . 8 . . 805 . . 35 . . 42 . 

X .485. .545, . 587. .625. . 653 , . 68 , . 7 . . 716 , .728, . 737 . . 745 . . 753 . . 754 . 

X .758, .76/ 

IDX=0 

DO  111  J=1 , 3 
DO  111  IL-1,17 
IDX=IDX+1 

IF(J-l)  111,109.110 

109  DUCT1 ( 1 , IL , 1 ) =TPX7 ( IDX ) 

DUCT1 ( 1 , IL . 4 ) =TP10 ( IDX) 

110  DUCT1 ( J . IL . 2 ) =TPX8 ( IDX) 

DUCT1 ( J . IL . 3 ) =TPX9 ( IDX ) 

DUCT1 ( J . IL . 5 ) =TP11 ( IDX ) 

111  CONTINUE 
C 

DATA  TP12  /2. .1.517,1.124.1.  .1. .1.  .1. .1.  .1. ,1. . 
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FILE  4 Continued 


X 1. .1. .1. ,1. .1. ,i. .1-/ 

DATA  TP13  /. 35. .565, .77. .81. .77. .71. .663. 

X .63. .606. .6, .61. .62. .635. .657. . 68 . . 7 , . 72 , . 35 , . 49 . . 685 . . 78 . . 76 . 

X .712.  .663, .63. .606. .6, .61. .62.  .635.  .657. . 68 . . 7 . . 72 . . 5 . . 61 , . 8 . 

X .825. .804. .74, .685. .657. .64. .64, .643, . 655 , . 67 . . 69 , . 715 . . 745 . .764/ 
DATA  TP14  /. 3. .56. .72, .845. .94,1.01,1.04, 

X 1.07.1.085.1.10.1.10.1.09,1.077,1.06.1.04.1.02.1.00, .3. .62. .836, 

X .91, .97.1.01,1.04.1.07,1.085,1.10,1.10.1.09,1.077.1.06.1.04.1.02. 
X 1.00, .767. .83, .907. .966.1.005.1.04.1.06,1.08.1.095,1.10.1,10, 

X 1.09.1.077.1.06.1.04,1.02.1.00/ 

DATA  TP15  /. 5, .5, .5. .5, .5, .5, .5. .5. .5. .5, .5, .5, 

X .5. .5. .5. .5, .5/ 

DATA  TP16  /. 65. .57, .54, .568. .645, .73. 

X .78. .78, .78. .78, .78. .78. .78. .78, .78. .78. .78. .485. .456. .468. .513. 

X .59. .655. .675. .675. .675. .675. . 675 , . 675 . . 675 . .675, .675. .675. .675, 

X .41, .41. .43. .476.  .536, .584, .607. .613. . 62 , . 62 . . 62 . . 62 , . 62 . . 62 . . 62 . 
X .62. .62/ 

IDX=0 

DO  116  J=l,3 
DO  116  IL=1 . 17 
IDX=IDX+1 

IF(J-l)  116.114,115 

114  DUCT2 ( 1 , IL, 1 ) = TP12 ( IDX) 

DUCT2 ( 1 , IL , 4 ) =TP15 ( IDX) 

115  DUCT2 ( J . IL . 2 ) =TP13 ( IDX ) 

DUCT 2 ( J . IL . 3 ) =TP14 ( IDX ) 

DUCT2 ( J, IL, 5 ) =TP16 ( IDX ) 

116  CONTINUE 
C 

DO  21  K=1 ,4,3 
DO  20  J=2 . 3 
DO  20  IL»1 , 17 

DUCTO ( J , IL . K ) = DUCTO ( 1 , IL , K ) 

DUCT1 ( J, IL. K) =DUCT1 ( 1 , IL, K) 

DUCT 2 ( J , IL . K ) = DUCT 2 ( 1 , IL . K ) 

20  CONTINUE 

21  CONTINUE 
C 

DATA  TP17  /-. 4, -.2,0. .0.2. 0.4. 0.6. 0.8.1./ 

DATA  TP18  /. 152. .169. .189. .213, .243. .290. .361, .450/ 

DATA  TP 19  /. 241. .248, .257, .266, .275. .277,  279, .293/ 

DATA  TP20  /. 41. .41. .41, .41, .41. .41, .41, .41/ 

DATA  TP21  /. 988. .983. .980,1.105. .955. .9487,1.039,1./ 

DATA  TP22  /. 5. .5. .5. .5. .5. .5. .5. .5/ 

DO  22  IL=1, 8 
SINVAL ( IL) =TP17 ( IL) 

BEND ( IL, 1 ) =TP18 ( IL) 

BEND ( IL , 2 ) =TP19 ( IL ) 

BEND ( IL . 3 ) =TP20 ( IL ) 

BEND ( IL . 4 ) =TP21 ( IL) 

BEND ( IL , 5 ) =TP22 ( IL) 

22  CONTINUE 
C 

DATA  TP23  /0 ., 0 ...  0215 ,. 08 ,. 1405 .. 1341 ,. 2673 .. 2646 .. 0944 
*, .3038, .443, .3922, .4468, .2464, .235. .2107, .2099, .1492. .1192. .1497. 

* .1202. .0957,0./ 

DATA  TP24  /0.,0.,0.,0.,0.,0.,0..0.,0.,0.,0..0..0.. 
♦0..0..0..0..0..0..0..0..0..1./ 

DATA  TP25  /I . 898E+6 , 2 . 204E+6 . 3 . 777E+6 , 7 . 885E+6 , 1 . 05E+7 . 

* 1.023E+7.1.027E+7.1.351E+7.5.54E+6.2.316E+7.4. 544E+7 , 1 . 054E+8 . 

* 2.215E+8,2.E+8,1.74E+8.1.549E+8,1.538E+8,1.071E+8,8.454E+7, 

* 1 . 009E+8 , 7 . 236E+7 , 6 . 108E+7 . 5 . 411E+7/ 

DATA  TP26  /2 . 326E+8 . 7 . 325E+7 . 8 . 596E+7 . 4 . 877E+7 . 5 . 477E+7 . 

* 4.258E+7,5.011E+7,6.085E+7,4.762E+7.6.542E+7.7.11E+7,7.11E+7, 

* 9,493E+7,1.813E+8,5.777E+8.8.534E+8.1.398E+8,6*0./ 

DATA  TP33  /0 . . 0 . , 0 . . 0 . , 0 . . 0 . . 0 . . 1 . 0 . 15*0 . / 

DO  26  M=1 , 23 
SEPI (M) =TP23 (M) 

STH (M) =TP24 (M) 

SINR (M) =TP25 (M) 

SGAM (M) =TP26 (M) 

SCOBLT(M) =TP33 (M) 

26  CONTINUE 
C 

DATA  TP27  /. 4359 ,. 324 .. 1566 .. 0722 ,. 0113 , 0 .. 1 . . 

* 3*0. , .1, .25. .25, .25. .15.0. ,1. ,3*0./ 

DATA  TP28  / . 9 . . 65 , . 4 . . 15 . 0 . . 1 . , . 9 . . 65 . . 4 , . 15/ 

IDX  = 0 

DO  28  1=1.4 
DO  28  J=1 , 5 
IDX=IDX+1 

AWTS ( I . J) =TP27 ( IDX) 

IF ( 1-2 ) 27,27.28 


FILE  4 Continued 


331 

332 

333 

334 

335 

336 

337 

338 

339 

340 

341 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352 

353 


27  AMOD ( I , J) =TP28 ( 1DX) 

28  CONTINUE 


DATA  TP29  /7 . E-9 . 7 . E-9 . 7 . 08E-9 . 6 . 72E-9 . 6 . 03E-9 . 

*5 . 43E-9 . 4 . 83E-9 . 4 . 48E-9 . 4 . 33E-9 . 4 . 23E-9 . 3 . 96E-9 . 3 . 3E-9 . 1 . 73E-9 . 

* 7.E-10,6.5E-10,6.07E-10,6. 72E-10 , 5.35E-10,3.88E-10,3.42E-10, 

* 3 . 27E-10 , 3 . 22E-10 , 3 . 2E-10/ 

DATA  TP30  /5.46E-9.5.13E-9.4.84E-9.4.61E-9.4.44E-9. 

* 4 . 13E-9 ,4.01E-9,3.39E-9,3.15E-9,3.09E-9,2.64E-9,1.97E-9,1.12E-9, 

* 0 . 229E-9 , 9*0 . / 

DATA  TP31  /22*0 . , 1 . / 

DATA  TP32  / 1 . 66E-9 , 1 . 46E-9 . 1 . 27E-9 , 1 . 08E-9 , 8 . 75E-10 , 

* 7 . 35E-10 , 6 . 44E-10 , 5 . 3E-10 , 4 . 45E-10 , 3 . 5E-10 , 2 . 56E-10 , 2 . 15E-10 , 

* 1 . 77E-10 , 1 . 22E-10 , 6 . 6E-11 , 3 . 9E-11 , 8 . 37E-11 , 6*0 . / 

DO  32  M=1 , 23 

SNYNEU(M) =TP29 (M) 

HNDRSN ( M ) =TP30 ( M ) 

THMLN(M) =TP31 (M) 

32  DGAM(M) =TP32(M) 

RETURN 

END 
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